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ABSTRACT 

A  new  data  processing  technique  is  suggested  to  estimate  the  de¬ 
lay  time  between  initially  uptraveling  energy  which  is  reflected  once 
at  the  earth  s  surface  and  initially  downtraveling  energy  on  earth¬ 
quake  seismograms.  The  method  uses  optimum  inverse  filter  together 
with  a  criterion  that  measures  the  simplicity  of  a  seismic  signal 
convolved  with  an  inverse  filter.  The  inverse  filters  are  designed 
to  extract  primary  energy  in  the  presence  of  surface  reflected 
energy  and  random  noise  on  the  basis  of  a  least-mean- square  error 
criterion.  Filter  design  is  dependent  on  the  delay  time  between 
primaries  and  surface  reflected  events,  their  amplitude  ratio  and  noise 
to  signal  power  ratio.  The  simplicity  criterion  was  devised  on  the 
assumption  that  a  maximum  in  the  concentration  of  normalized  seismic 
signal  energy  above  a  minimum  level  indicates  which  filter  was  most 
correctly  designed.  This  is  visualized  as  an  expression  of  the  hypo¬ 
thesis  that  the  primary  seismogram  is  generated  by  a  few  large  dis¬ 
continuities,  rather  than  by  many  minor  boundaries. 

The  technique  was  applied  to  band-limited  synthetic  signals  that 
contained  several  primary-secondary  pairs  in  the  presence  of  random 
noise.  Of  27  synthetic  signals  which  were  analyzed,  the  procedure 
successfully  selected  the  correct  delay  time  in  22  cases, 

Four  actual  earthquake  seismograms  were  then  analyzed.  The 
procedure  selected  a  delay  Lima  for  each  earthquake.  Focal  depths 
computed  from  the  selected  delay  times  appeared  quite  reasonable  when 


lil 


compared  with  depths  for  the  same  earthquakes  published  by  the  United 
States  Coast  and  Geodetic  Survey.  In  all  cases,  the  inverse  filter  de¬ 
signed  for  the  selected  delay  time  considerably  simplified  the  original 
seismogram. 

It  is  concluded  that  the  technique  provides  a  reasonable  estimate 
of  the  delay  time  between  primary  and  surface  reflected  energy. 
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INTRODUCTION 

General  Statement  of  the  Problctn 

This  thesis  is  crncerned  with  the  problem  of  extracting  infor¬ 
mation  from  a  bandlimited  signal  in  the  presence  of  corrupting  noise. 
In  particular  the  signal  is  a  seismic  signal  generated  by  an  earth¬ 
quake  at  or  near  the  earth’s  surface.  The  information  desired  is  a 
correct  estimate  of  the  focal  depth  of  the  earthquake's  source. 

An  accu’^ate  estimate  of  the  source  depth  would  help  to  differ¬ 
entiate  natural  earthquakes  from  clandestine  nuclear  explosions,  .since 
the  depth  of  a  nuclear  explosion  is  .subject  to  practical  limitations. 
For  example,  all  disturbances  with  focal  depths  greater  than  say  10  km 
could  be  classified  as  natural  with  reasonable  certainty. 

The  Coast  and  Geodetic  Survey  is  the  agency  in  the  United  States 
which  routinely  collects  seismic  data  from  stations  throughout  the 
world.  They  compute  the  focal  depth  of  an  earthquake  in  kilometers 
below  the  mean  sphere  based  on  P-arrivals  and  on  the  Jef freys-Bul len 
travel-time  tables  of  1958.  The  survey  made  the  following  statement 
concerning  the  accuracy  of  their  published  focal  depths. 

"A  freely  determined  depth  is  checked  against  depth 
phases  interpreted  from  available  seismograms  or 
reported  by  cooperating  observatories.  Depths  may 
be  restricted  to  agree  with  these  depth  phases  if 
agreement  tc  within  the  stated  accuracy  is  not 
obtained.  In  the  case  of  shallow  earthquakes  or 
smaller  earthquakes,  the  exact  depth  cannot  usually 
be  determined  precisely  and  the  depths  quoted 
represent  a  judgnent."  (U.  S,  Department  of 
Commerce  Coast  and  Geodetic  Survey,  1964) 
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A  fo^al  depth  tar  an  earthquake  published  by  thi 


s  survey  may  be  in 


error  bv  as  much  t-j ^ morose  j  t-  ,  , 

,s  ^^^^umeters.  lUinst  and  Engdahl,  1962). 


Frevious  Studies 


several  investigators  have  tried  with  varying  degrees  of  success 
to  develop  techniques  that  will  obtain  more  accurate  focal  depths  from 
seismograms.  Most  techniques  or  suggested  methods  fall  in  one  of  three 
categories:  visual  recognition  of  seismic  events  with  arrival  times 
related  to  focal  depth,  spectral  analyses  and  linear  filter  operations. 

Visual  recognition  is  the  most  obvious  technique  and  simplest  to 
apply.  In  1936  Gutenberg  and  Richter  suggested  using  visible  and 
clearly  identifiable  pulses  such  as  P,  5  and  pP  (Howell,  1959).  The 
paths  followed  by  these  pulses  are  illustrated  in  Figure  1.  S-waves 
follow  essentially  the  same  path  as  P-waves  but  at  a  slower  velocity. 

If  the  velocity  profile  at  the  source  is  known,  the  delay  time 
(difference  in  arrival  times)  between  pP  and  P  gives  a  direct  estimate 
of  the  focal  depth.  A  knowledge  of  the  epicentral  distance  and  time  of 
origin  of  the  event,  combined  with  arrival  times  for  P  and  S,  also 
gives  a  direct  estimate  of  the  depth  of  focus.  For  shallow  earthquakes 
sP  may  be  easier  to  identify  than  pP.  This  is  due  to  the  fact  that  sP 
travels  to  the  surface  at  a  slower  velocity  than  pP,  and  hence  is 
separated  by  a  greater  time  difference  from  P.  Kondorskaya  (1956)  has 

in  some  cases  identified  sP  and  used  its  arrival  time  to  determine  the 
focal  depth  of  a  shallow  event. 
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EPICEfITRAL  DISTANCE 


RECEIVER 


Fig.  1.  Gutenberg  and  Richter  ray-path  configuration. 


For  crustal  earthquakes  other  pulses  may  be  used  provided  that 
the  receiving  station  is  at  a  distance  (D)  from  the  source  that  lies 
within  the  first  zone.  This  zone  includes  distances  up  to  1000  kilo¬ 
meters,  for  which  the  curvature  of  the  earth  can  be  neglected.  When 
this  is  the  case,  (i.e.,  D  is  known  and  D<1000  km)  Thirlaway  (1961) 
suggests  using  the  delay  time  between  and  or  and  to 
determine  the  focal  depth.  The  ray-paths  for  these  pulses  are  illus¬ 
trated  in  Figure  2.  De  Bremaeker  (1955)  has  shown  that  the  variation 

of  P  amplitude  with  epicentral  distance  is  also  related  to  the  source 
n 

depth  of  crustal  earthquakes. 


< - D<1000  km  - — — -*■  RECEIVER 


2,  Thirlaway's  ray-path  configuration  for  crustal 
earthquakes  recorded  in  the  first  zone. 


. . .  . .  . . .  . .  .  ■!  ,1.  ^ 
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All  these  techniques  depend  on  visual  identification  of  clearly 
recognizable  pulses.  This  is  a  severe  restriction  when  focal  depths 
are  shallow  and  consequently  differences  in  arrival  time  are  so  small 
that  pulses  overlap.  Furthermore  background  noise  level  impairs  visual 
identification. 

As  an  alternative  to  direct  visual  recognition  of  seismic  events, 
correlation  procedures  have  also  been  suggested.  Time  autocorrelation 
of  the  seismic  signal  was  used  by  The  California  Research  Corporation 
(Cal.  Research)  (1961)  and  Rothman  (1964)  to  estimate  the  delay  time 
between  two  sets  of  pulses  whose  time  separation  is  related  to  the 
source  depth.  This  process  measures  the  correlation  of  the  signal  as 
a  function  of  time  as  it  is  shifted  across  itself.  Mathematically 
this  is  expressed  as 

/■*oo 

«P(T)  =  1  x(t-T)x(t)dt 

J~(30 

where  q5(T)  is  the  time-autocorrelation  function,  and  T  is  the  time 
displacement  (shift)  of  x(t)  (Lee,  1963).  The  maximum  correlation 
occurs  when  T  =  0  and  is  positive.  For  a  seismic  signal  containing  two 
pulses  or  two  sets  of  pulses  having  the  same  shape  but  differing  in 
amplitude,  there  will  be  a  relative  maximum  in  the  absolute  value  of 
cp(T)  at  T  equal  to  their  time  separation.  If  the  pulse  sets  have  the 
same  polarity  this  correlation  will  be  positive^  and  if  they  tiave 
opposite  polarities  it  will  be  negative.  When  more  than  two  pulse 
sets  are  included  in  the  analysis,  other  correlation  peaks  or  troughs 
will  appear  in  the  function.  If  this  is  the  case,  the  position  of  the 
relative  maximum  may  or  may  not  correspond  to  the  correct  delay  time. 
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The  difference  frequency  between  minima  in  the  amplitude  spectrum 
of  a  transient  containing  two  pulses  or  two  functions  of  the  same 
shape,  has  been  shown  by  Nakam.ura  (1963)  to  be  related  to  their  delay 
time.  When  the  pulses  are  out  of  phase  Nakamura  showed  that  minima 
in  the  spectrum  will  appear  at  integer  multiples  of  the  harmonic  fre¬ 
quency  given  by 

^o=  I 

where  T  is  the  time  separation  of  the  two  pulses. 

Rotliman  (1964)  obtained  spectra  for  two-dimensional  model  re¬ 
cords,  and  was  successful  in  picking  the  difference  frequency  related 
to  the  focal  depth  when  the  duration  of  analysis  included  only  the  two 
pulses  of  interest.  Some  ambiguity  in  picking  the  correct  difference 
frequency  resulted  when  more  pulses  were  included  in  the  analysis. 

That  is,  the  difference  frequency  does  not  uniquely  deteririine  the  de¬ 
lay  between  two  pulses  unless  they  are  tiie  only  major  pulses  in  the 
signal . 

Bogert  et  al  (1963)  used  the  time  autocorrelation  function  of  a 
seismic  signal  to  find  the  Cepstrum.  Ihey  defined  the  Ccpstrum  to  be 
the  spectrum  of  the  log-power  spectrum.  By  band-broadening  the  log- 
power  .spectrum  of  the  autocorrelation  function  they  found  that  they 
could  sometimes  locate  periodicities  (related  to  the  delay  time)  in  the 
Cepstrum  that  were  not  evident  in  the  autocorrelation  function.  Their 
results  were  quite  dependent  on  the  type  of  high-pass  filter  they  used 
for  band-broadening  the  log-power  spectrum. 
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Both  Keylis  Borok  (1962)  and  Cal.  Research  (1961)  used  the 
amplitude  spectrum  of  surface  waves  to  differentiate  between  earth¬ 
quakes  and  nuclear  explosions.  Energy  distribution  as  a  function  of 
frequency  actually  depends  on  several  factors,  one  of  which  is  the 
source  depth.  However,  more  predominant  factors  such  as  detector  lo¬ 
cation,  earth  layering  and  source-detector  separation  are  more  im¬ 
portant  in  determining  the  frequency  distribution  of  energy,  thus 
limiting  the  power  of  this  technique. 

To  be  successful  most  of  these  procedures  can  only  be  applied  to 
relatively  broad-band  signals.  However,  seismic  signals  are  inherently 
narrow-band  (2-3  octaves).  This  is  because  the  amount  of  attenuation 
of  seismic  energy  due  to  absorption  (earth  filtering)  is  frequency  de¬ 
pendent.  The  result  of  this  is  a  rapid  loss  of  high-frequency  energy 
in  seismic  waves  with  distance  from  the  source  (Howell,  1959).  in 
addition,  the  instrumentation  at  the  receiver  and  the  local  geology  at 

both  the  source  and  receiver  tend  to  limit  the  bandwidth  of  seismic 
energy. 

Recognition  of  th''  facts  has  stimulated  researchers  to  attempt 
to  band-broaden  seismic  data.  Sufficient  band -broadening  should 
visually  simplify  a  seismic  signal  and  help  to  make  individual  pulses 
recognizable.  At  the  Massachusetts  Institute  of  Technology  under 
the  direction  of  Simpson  (1961,  1962)  and  Cal.  Research  (1961),  re¬ 
search  has  been  directed  towards  finding  an  estimate  of  the  earth 
filter.  Convolution  of  the  inverse  of  the  earth  filter  with  the 
seismic  signal  should  yield  a  broader-band  signal.  Scientific  reports 
from  M.I.T.  (Simpson,  1961,1962)  and  Cal.  Research  (1961)  suggest  that 


tne  sha.Tt-  n;  t’le  first  primary  paise  in  a  snismic  signal  would  be  a 
reasui."'bl\  good  estimate  of  the  earth  filter.  However,  seismograms 
filtered  with  tiie  invetEe  -f  t;ie  primary  pulse  were  not  significantly 
simplified.  This  inability  to  simplify  the  seismogram  with  inverse 
filtering  of  the  earth  filter  is  due  to  several  reasons.  In  the  first 
place,  the  process  for  finding  the  inverse  filter  is  extremely 
sensitive  to  truncation  or  smooching  of  the  original  pulse  shape 
(Watson,  196h  personal  communication)  Furthermore,  it  is  very  diffi¬ 
cult  to  obtain  a  good  estimate  of  the  earth  filter  from  the  seismic 
signal,  due  to  corrupting  noise  and  overlapping  pulses.  Even  if  the 
earth  filter  were  known  reliably  it  is  doubtful  that  much  band¬ 
broadening  could  be  achieved,  due  to  the  presence  of  noise  and  the 
limited  bandwidth  of  the  seismic  filter  (Howell  et  al,  1953). 

Due  to  the  problems  innerent  in  the  application  of  many  of  tne 
techniques  just  discussed,  it  is  apparent  that  a  new  apnroach  to  the 
focal  depth  problem  is  necessary  In  the  following  pages  we  will  in¬ 
vestigate  a  new  technique  chat  circumvents  many  of  the  difficulties 
discussed  in  tnis  section. 

Explanation  of  the  Suggested  Procedure 

let  us  consider  a  seismic  disturbance  at  some  arbitrary  depth 
beneatn  tlie  earth's  surface.  The  disturbance  emanates  energy  in  ail 
directions.  Energy  sent  initially  downward  travels  directly,  as  well 
as  through  a  combination  cf  refractions  and  reflections,  to  a  receiver. 
At  the  receiver  tiiis  energy  creates  a  series  of  pulses  on  the  seismo¬ 


gram  wtiich  we  call  primary  events. 


Energy  emanated  initially  upward  is 


8 

reflected  at  the  earth's  surface  (or  from  intermediate  discontinuities), 
and  then  follows  travel  paths  to  the  receiver  similar  to  those  taken  by 
the  initially  downward-t ravel ing  energy.  This  energy  also  creates  a 
series  of  pulses  on  the  seismogram  which  we  call  secondary  events  or 
ghosts  In  addition,  other  events  such  as  surface  waves  and  micro¬ 
seisms  are  recorded  These  events  and  all  others  which  are  neither 
primaries  nor  ghosts  constitute  noise  in  this  mathematical  model 
(Howell  et  al ,  1963) . 

Let  us  suppose  chat  eacn  primary  event  has  an  associated  ghost. 

This  ghost  is  presumed  to  have  the  same  pulse  shape  as  the  rimary  but 

IS  reversed  in  polarity.  The  ratio  of  ghost  to  primary  amplitude  (R  ) 

o 

IS  determined  primarily  by  the  surface  reflection  coefficient,  and 
attenuation  of  ghost  energy  in  its  longer  travel  path.  Due  to  its 
longer  travel  path,  the  ghost  lags  behind  the  primary  by  a  time  1  de¬ 
pendent  on  the  depth  of  the  disturbance. 

Now  if  we  could  design  an  operator  that  could  suppress  the 
secondary  events  on  a  seismogram  (inverse  filtering;,  only  the  primary 
events  would  remain  in  the  presence  of  filtered  noise  Presumably  the 
best  operator  would  be  designed  with  an  exact  knowledge  of  T  and  R 

o ' 

fc  be  effective  it  would  also  have  to  operate  in  the  presence  of 
corrupting  noise.  However,  T  and  in  general  would  not  be  known.  If 
T  were  known  there  would  be  no  problem.  But  if  several  inverse  oper¬ 
ators  were  designed  covering  the  range  of  possible  values  of  R  and  I, 

o 

one  of  these  would  be  most  nearly  correct,  at  least  to  within  the 
intervals  of  variation  chosen  for  and  T.  To  ascertain  which  inverse 

Operator  and  filter  will  be  used  as  synonyms  in  this  paper. 
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operator  based  on  particular  values  of  and  I  was  correct,  it  would 
be  necessary  to  determine  wnlch  Inverse  filtered  seismogram  most  nearly 
contains  only  primary  events.  A  criterion  wPicti  attempts  to  measure 
the  simplicity  of  the  inverse  filtered  seismogram  could  be  devised  to 
make  the  latter  determination.  The  primary-ghost  separation  would  tf  - 
be  known  and  hence  the  depth  of  the  source  could  be  obtained,  if  the 
velocity  profile  were  known. 

The  procedure  just  outlined  would  not  require  a  knowledge  of  the 
earth  filter.  Further,  it  should  be  reasonably  independent  of  the 
number  of  primary-ghost  pairs  provided  that  the  pulse  shape  of  each 
ghost  is  the  same  as  its  associated  primary,  and  that  their  amplitude 
ratio  is  the  same  for  all  pairs  included  in  the  analysis. 

Mathematical  Model  for  the  Present  Study 

Mathematically  the  above  ideas  may  be  expressed  as  follows  Let 
the  pulse  shape  of  the  events  that  initially  traveled  downward  be  b(t). 
In  this  model,  b(t)  includes  the  shot  generation  pulse,  the  earth 
transmission  attenuation  characteristic  and  the  filtering  due  to  re¬ 
cording  instrumentation.  It  is  equivalent  to  the  b(t)  function  as  de¬ 
fined  by  Sengbush  ec  al  (1961)  with  additional  filtering  due  to 
attenuation.  By  writing  b(t)  instead  of  b(t,T),  where  T  is  travel 
time,  we  have  assumed  time-invariance  of  the  seismic  wavelet,  and  thus 
have  asserted  that  the  attenuation  is  identical  for  all  events  on  a 
given  sei..j30gram.  By  defining  b(t)  in  this  way  the  primary  events  on 
the  seismogram  may  be  expressed  as, 


P(t)  b(t)*r(t) 


(1) 
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(Simpson,  1961),  where  r(L)  is  the  response  of  the  earth  to  a  plane-wave 
impulsive  source  propagating  without  attenuation  and  recorded  with  per¬ 
fect  instruments.  r(t)  includes  the  direct  pulses  (for  teleseismic 
events),  reflections  and  refractions  which  have  ghost  counterparts. 

All  other  events  not  having  ghost  counterparts  such  as  surface  waves 
are  considered  part  of  the  noise.  This  r(t'  is  analogous  to  the 
"reflectivity  function"  introduced  by  Peterson  et  al  (1955)  which  is 
restricted  in  their  model  to  include  only  reflected  events.  The 
asterisk  .stands  for  convolution 

x(c)*y(t)  =  \  x(T)y(t-T)dT.  (2) 

I f  we  consider  only  those  ghosts  which  reflected  from  the 
earth's  surface  (these  are  believed  to  contain  the  most  energy),  then 
the  secondary  events  may  be  exprec.‘’ed  by 

G(t)  =  -R^b(t-T)*  r(t).  (3) 

riere  is  the  ratio  of  ghost  to  primary  amplitude  nd  may  range  from 
almost  zero  to  slightly  greater  than  unity  in  magnitude.  I  is  the 
time  required  for  energy  to  travel  to  the  earth's  surface  from  the 
source,  and  back  down  again  to  the  same  depth  as  the  disturbance. 

The  expression  for  tl.e  recorded  signal  is  the  sum  of  the  primary 
(equation  (1)),  and  secondary  (equation  (3))  events  together  with 
noise. 


S(t)  =  P(t)+G(t)+n(t). 


(4) 
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Noise,  n(f.),  constitutes  all  other  e\^ents  which  are  neither  primaries 
nor  ghosts.  Rewriting  (4)  we  have, 

S(t)  =  b(t)*  r(t)-R  b(t-T)*r  (t)-n(t).  (5) 

o 

This  may  be  written, 

S(t)  =  b(t)^  r(t)^  [<^(t)-R^c^  (t-T)j -n(t).  (6) 

^(t)  is  the  impulse  function 

{^(t)  =0  t  0 

($(c)dt  =  1  t  ^  0 

-  00 

(Churchill,  1958).  (^(t)  has  tne  additional  property  that 

S  (t)*y(c)  ,=  y(t) . 

Now  let 

h(t)  =  ^(t)  -  R^5(t-r)  .  (7) 

h(t)  is  a  doublet  function  consisting  of  two  impulses,  one  positive  and 

the  other  negative,  and  separated  in  time  bv  an  amount  I.  R  is  the 

o 

ratio  of  their  amplitudes.  Ibis  doublet  function  will  be  called  the 
ghost  filter.  Subs  'tuting  for  h(t)  in  (6),  we  have, 

S(c)  b(t)i^  r(t)>-h(t)  *  n(t)  (8) 

S(t)  in  (8)  is  the  proposed  mathematical  expression  for  the  recorded 
sequence  of  events  generated  by  a  source  at  an  arbitrary  depth  beneath 
the  earth's  surface  Figure  3  is  a  schematic  of  ttie  proposed  data 


r 
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generation  process. 


S(t)  =  F(t)^G(t)+n(t) 


Fig. 3.  Schematic  of  the  proposed  data  generation 
hypothesis . 


In  summary  the  data  generation  hypothesis  is  based  on  the 
following  assumptions.  All  primary  events  generated  by  a  seismic 
disturbance  pass  through  a  ghost  filter.  This  filter  introduces  a 
pulse  (ghost)  after  each  primary  that  has  the  same  pulse  stiape  as  the 
primary  but  is  reversed  in  polarity.  It  is  delayed  by  a  time  that  is 
related  to  the  depth  of  the  disturbance.  The  amplitude  ratio  of  the 
primary  to  ghost  is  fixed  for  all  primary-ghost  pairs  on  a  given 
seismogram.  All  ocher  events  which  are  neither  primaries  nor  ghosts 
constitute  noise. 


Focal  Depth  Estimation  Procedure 

Suppose  and  T  are  known,  then  h(t)  is  completely  specified. 

Now  devise  an  operator  h  ^(t)  which  has  the  property 

h‘^(t)^h(t)  Si  Sit).  (9) 

This  is  an  approximation  since  h  ^(t)  must  be  designed  to  extract  the 
primary  events  from  a  signal  in  the  presence  of  noise,  and  hence  a  per¬ 
fect  inverse  operator  or  exact  inverse  is  not  required  (Foster  et  al, 
1962).  It  follows  that, 


(10) 


h  ^(t)*S(t)  h  ^(t)*b(t)*  r(t)*h(t.)  +  n'(t) 


>ss  b(t)^  r(t)  -  n'(t) 


P(t)  +  n'(t) 


where  n'(t)  h  (t)*n(t).  h  (t)  is  the  mathemaLical  inverse  of  h(t), 
a  doublet  function  in  the  sense  to  be  described  under  Theoretical  De¬ 
sign  of  Inverse  Operators.  The  design  of  h  ^(t)  requires  a  knowledge 

of  R  and  T.  Since  in  general,  knowledge  of  these  parameters  will  be 
o 

limited,  it  is  suggested  that  an  array  of  inverse  operators  h  (t,R^j, 

I  )  be  generated  for  various  values  of  R  and  T.  The  j  and  k  sub- 
R  O 

scripts  on  R  .  and  T,  refer  to  possible  values  of  the  R  and  T  parame- 
^  oj  k  o 

ters . 

The  given  earthquake  seismogram  is  filtered  by  each  of  these  in¬ 
verse  operators  (Figure  4)  producing  a  suite  of  inverse  filtered 

traces  P.,  ,  each  for  a  different  R  . ,  T,  combination  in  tne  inverse 
jk  oj  k 

filter.  That  combination  of  R  .  and  T,  for  the  inverse  filter  which 

OJ  k 

yields  the  "best”  estimate  of  the  primaries  is  considereo  the  most 

correct  combination  for  the  signal  analyzed.  By  "best"  we  mean  that 

set  of  primaries  whose  energy  is  concentrated  in  the  fewest  number  of 

large  pulses.  This  set  of  primaries  relative  to  the  other  sets  should 

have  the  simplest  structure.  The  simplicity  criterion  is  an  attempt  to 

quantify  this  concept.  This  criterion  will  yield  an  estimate  of  I  and 

R  ,  called  T  and  R  ,  T  together  with  an  appropriate  velocity  profile 
o  o 

will  be  used  to  estimate  the  deptn  of  focus  of  the  earthquake.  The 
complete  estimation  procedure  is  illustrated  schemal ical iy  in  Figure  4. 
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P(t:;T,R  ) 

O 


Pig-  4.  Schematic  of  the  proposed  focal  depth  estimation  procedure. 

Theoretical  Design  of  the  Inverse  Operators:  An  expression  for 
the  exagi  inverse  operator  for  the  doublet  function  h(t)  was  derived  by 
Lindsey  (1960),  and  has  the  form, 

00 

h'^(t;R^,T)  =  (t-nT). 

n=o 

This  expression  for  h  ^(t)  is  neither  practical  nor  desirable.  First 
of  all,  it  requires  an  infinite  number  of  points  and  secondly,  it  is 
inadequate  for  a  system  containing  noise.  The  inverse  filter  we  re¬ 
quire  should  be  capable  of  doing  an  effective  job  in  the  presence  of 
noise.  It  should  also  contain  only  a  finite  number  of  points  so  that 
it  may  be  used  in  a  computer. 

The  general  technique  for  designing  an  optimum  filter  was  de¬ 


veloped  by  Norbert  Weiner  (1949).  He  formulated 


a  criterion  or  least- 


square  measure  that  gave  an  objective  ir.dication  of  filter  performance. 
An  optimum  filter  is  a  filter  whicn  minimizes  the  mean-square  error 
difference  between  a  desired  output  and  the  actual  output  of  the  filter 
as  explained  below.  In  addition,  Weiner  statistically  described  the 
signal  and  noise. 

The  procedure  for  finding  the  optimum  inverse  f 11  ter  was  outlined 
by  Foster  et  al  (1962).  The  problem  was  formulated  in  the  manner  indi¬ 
cated  by  the  flow  chart  in  Figure  5. 


Fig.  5.  Flow  chart  for  formulation  of  the  optimum  inverse  filter 
problem. 

The  purpose  was  to  determine  c(t)  which  is  in  a  certain  sense  the  in¬ 
verse  to  the  filter  labeled  a(t).  c(t)  is  restricted  to  be  a  linear, 
band-limited,  time- invariant,  finite  memory  filter.  The  data  that  c(t) 
operates  upon  is  the  filtered  signal  mixed  with  noise 

i+aO 

a(t)s(t-I)dT  ^  n(t).  (11) 

'00 

The  output  of  c(t)  is  an  estimate  of  s(t)  and  is  given  by 

( 

^  \ 

d(t)  \ 


c( t) i( t-T)dT, 


(12) 
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where  "  is  used  to  indicate  that  d{t) 
error  made  in  estimating  s(t)  is 


is  an  estimate  of  s(t). 


The 


e(t)  -  s(t)  -  d(t)= 

According  to  Weiner  (1949)  the  best  estimate  d( t)  of  s(t)  is  obtained 
when  a  filter  c(t)  is  used  that  minimizes  the  mean-square  error.  The 
mean-square  error  in  this  case  is  given  by 


lim  _1 
T-**oo  2T 


(14) 


For  the  case  when  only  a  minimum  knowledge  about  the  statistical 
properties  of  the  signal  and  noise  is  available  (maximum  entropy  case) 
Foster  made  the  following  assumptions: 

1.  The  signal  s(t)  and  noise  n(t)  are  uncorrelated 
random  processes. 

2.  The  autocorrelation  function  of  s(t)  is 


and  for  n(t) 

c/  s  (I) . 


(15) 


i  2 

^s  signal  power  and  is  the  noise  power. 

2  2 

The  ratio  /G^  is  called  the  noise/signal  power 


ratio  (R). 


Using  these  assumptions  Foster  derives  a  stable  optiiaum  inverse  filter 
e(t)  that  minimizes  the  Weiner  mean-square  error.  eft)  was  chosen  irom 
the  class  of  all  band-limited,  linite  memory  functions. 

Using  Foster's  procedure  an  optimum  inverse  filter  can  be  derived 
for  the  doublet  function  h(t)  for  any  combination  of  For  this 

case  the  flow  chart  of  the  problem  is  illustrated  in  Figure  6.  Here 
we  desire  an  optimum  inverse  filter  h  ^(t)  that  will  minimize  the  mean- 
square  error 


1  im  1 
T->oo  2T 


(16) 


where  z(t)  is  an  estimate  of  the  impulse  function  d(t). 


Fig.  6,  Flow  chart  of  problem  for  finding  an  optimum  inverse  filter 
for  the  doublet  function. 


In  using  the  Foster  model  for  derivation  of  the  inverse  filters, 
we  have  assumed  that  the  seismic  noise  is  wide  sense  stationary  and 
white.  These  assumptions  are  probably  not  true  fot  the  seismic  data 
generation  process,  but  are  nevertheless  appropriate  when  no  further 
knowledge  concerning  the  noise  statistics  are  available  (Foster  et  al. 


1964) . 
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A  typical  inverse  operator  together  witn  the  input  and  output 
signal  is  shown  in  Figure  7.  ine  interval  between  sample  points  is 
equal  to  the  assumed  delay  time  between  the  primary  ant  ghost.  Con¬ 
volution  in  the  time  domain  with  the  doublet  results  in  a  single 
positive  impulse  together  with  several  small  negative  impulses 
(secondary  spikes).  These  negative  impulses  are  generated  since  the 
inverse  operators  are  only  designed  to  do  the  best  job  in  the  presence 
of  noise,  under  t^.e  finite  memory  and  band-limit  restrictions. 


h(t)*b'^(t)  =  Z(t) 

Fig.  7.  Output  signal  z(t)  found  by  convolution  of  h  \t) 

(correct  inverse)  with  the  input  signal  h(t). 

Physically  the  convolution  of  h( t)  with  h'^(t)  proceeds  as 
follows;  The  initial  impulse  in  h(t)  reproduces  h"^(t).  The  second 
impulse  in  h(t),  (i.e.,  -R^,$(t-T))  reproduces  h  ^(t)  multiplied  by 
'^o  delayed  T  seconds.  Thus  the  second  impulse  produces 
-R^h  (t-T).  Finally,  h  ^(t)  and  -R^h  ^(t-T)  are  summed  to  produce  the 
function  labeled  h(t)  *  h’^(t)  in  Figure  7. 

In  the  frequency  domain  convolution  corresponds  to  multiplication 
of  spectra.  One  period  of  the  amplitude  spectrum  of  the  ghost  filter 


h(t) 


h‘^(t) 


Iv 

h(t)  and  the  amplitude  spectrum  of  the  appropriate  inverse  h  ^(t)  are 
shown  in  Figure  8a.  The  product  of  these  spectra  yields  a  spectrum 
that  is  broadei -band  than  the  ghost  filter,  i.e.,  more  nearly 
approaches  the  spectrum  of  acfCt)  function). 

To  examine  further  the  effects  of  inverse  filtering  on  a  signal 
in  the  frequency  domain,  a  signal  containing  three  positive  impulses 
(separated  in  time)  was  constructed  and  subsequently  filtered  to 
simulate  real  data. 

Figure  8b  is  the  amplitude  spectrum  of  three  positive  primary 
impulses  that  have  not  been  filtered.  Below  this  is  the  amplitude 
spectrum  of  the  primaries  filtered  by  the  ghost  filter  and  b(t).  The 
next  amplitude  spectrum  in  Figure  8c  was  obtained  using  the  correct 
inverse  filter  for  the  giiost  filter.  The  final  spectrum  was  obtained 
using  an  inverse  filter  designed  for  a  shorter  primary-gliost  delay 
time.  It  is  quite  clear  that  this  filter  boosted  both  the  low  and  high 
frequencies  in  the  spectrum.  That  is,  the  deconvolved  signal  using  the 
wrong  filter  is  broader-band  chan  the  signal  produced  by  t.he  correct 
inverse  filter.  The  correct  inverse  filter  did  band-broaden  the  signal 
but  not  as  much  as  the  incorrect  filter.  In  a  few  cases,  this  over¬ 
compensation  phenomenon  is  thought  to  be  one  reason  for  the  failure  of 
the  simplicity  criterion  to  differentiate  between  two  sets  oi  estimated 
primaries  P.  This  criterion  will  be  explained  more  fully  in  the  next 


section . 


20 


FREQUENCY  IN  CPS 

Fig.  8a.  One  period  of  The  amplitude  spectra  of  the  ghost 
filter,  ghost  Inverse,  and  their  product. 


H 
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Fig.  8b. 


Amplitude  spectra  of  input  signal  before  inverse 
filtering. 


Fig.  8c.  Amplitude  spectra  of  input  signal  after  inverse  filtering 
(h^^(f)  correct  inverse,  h”^(f)  incorrect  inverse).  The 
variable  f  is  used  to  indicate  that  the  functions  above 
are  Fourier  transforms  of  time  functions,  e.g.,  the 
Fourier  transform  of  r(t)  is  designated  by  r(f). 


Pritnary-es timation  Method;  Ideally,  we  would  like  co  choose  one 
of  the  inverse  filtered  records  as  being  the  best  estimate  of  the  pri- 

A 

mary  events  P(t).  As  explained  on  page  13  ,  this  would  yield  T  and  R 

o, 

estimates  of  T  and  R  respectively.  How*ever,  we  do  not  know  the  actual 

o 

primary  function.  Therefore  we  must  use  a  criterion  independent  of  the 
true  set  of  primaries  to  make  this  choice.  We  assume  that  the  inverse 
filtered  record  (contained  in  a  suite  cf  inverse  filter\_d  records) 
which  contained  the  most  nearly  correct  set  of  primaries  would  be  the 
simplest  in  visual  character.  Our  reasoning  here  comes  from  the  re¬ 
sult  of  convolution  of  the  actual  h(t)  filter  present  in  a  seismogram, 

with  an  h  ^(t)  filter  based  on  incorrect  values  for  R  and  T.  Figure  9 

o 

show  one  such  possible  incorrect  h  \t)  convolved  with  h(t).  Compari¬ 
son  of  the  output  of  this  operation  with  the  desired  operation  shown  in 


h(  t) 


h’^(t) 


z(t) 


Fig.  9.  Output  signal  z(t)  found  by  convolution  of  h  (t) 
(incorrect  inverse)  with  input  signal  h(t). 


Figure  7  indicates  that  the  correct  inverse  produces  a  single  large 
impulse  (with  attendent  small  negative  spikes)  whereas  the  incorrect 
inverse  produces  several  large  impulses  distributed  over  a  longer 


time  interval.  Thus  if  we  can  measure  the  concentration  ot  energy  per 


unit  time,  i.e. ,  quantify  the  visual  difference  between  Figures  7  and 
9,  we  would  have  an  objective  criterion  for  choosing  the  result  of 
Figure  7  over  chat  in  Figure  9. 

As  a  measure  of  energy  concentration  per  unit  time  we  choose  the 
function 


Ciu;R^j.Tk) 


i  .  1,2,. . .n  (17) 


u  $  P(t;R  .,T.  ) 

OJ  K 


given  in  equation  (17).  Here  P(t;R  ,»T.)  is  the  deconvolved  seismo- 

Oj  K 

gram  based  on  j »  ^  arbitrary  but  fixed  level,  A^  is 

A 

the  total  time  for  which  the  magnitude  of  P  is  greater  than  or  equal 

/V. 

to  u.  The  lower  integral  represents  the  total  energy  in  P  for  u  -  - . 
Figure  10  illustrates  how  C  is  determined-  The  C(u;R  .»T,  )  function 

O  j  K 

is  defined  for  various  amplitude  levels  u,  and  thus  represents  the 

A 

energy  of  P  above  the  value  u  per  unit  time,  i.e.,  the  energy  concen¬ 
tration  above  a  given  u. 
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Fig.  10.  Diagram  of  the  parameters  involved  in  the  formulation  of 
the  simplicity  criterion. 


We  choose  as  the  correct  estiiuate  of  P(t),  that  P{t,R  ,T  )  wliich 

oj  k' 

maximizes  G(u;R  ,,T,  )  for  the  smallest  value  of  u,  C(u.R  ,T,  )  as  a 

Oj  K.  OJ  K 

function  of  u  will  have  a  maximiiir.  for  a  single  inverse  filtered  pri¬ 
mary-ghost  impulse  pair  when  u  just  exceeds  the  magnitude  of  the 
largest  secondary  spike  (for  example  the  output  h  ^(t>-eh(t)  of  Figure 
7).  Thus  C(u;R  . ,T,  )  will  decrease  bevond  this  value  of  u.  It  will 

OJ  k 

also  be  clear  from  Figure  9  that  a  maximum  in  concentration  for  z(t) 

will  occur  for  a  larger  value  of  u,  and  furthermore  the  maximum  will 

not  be  as  great  as  that  for  z(t)  in  Figure  7.  This  is  because  the  A. 

in  the  denominator  of  C(u.;R  .,T.  )  for  z(t)  in  Figure  9  will  be  large 

o  j  k 

compared  to  the  for  z(t)  in  Figure  7.  Appendix  A  presents  a  proof 
that  a  choice  of  primaries  based  on  this  criterion  is  va^id  wnen  the 
seismogram  consists  of  a  single  impulse  without  noise,  and  the  inverse 
filter  procedure  is  exact. 

Part  of  the  ooject  of  this  thesis  is  Co  test  the  validity  of 
this  same  criterion  when  applied  to  more  realistic  data  such  as 
synthetic  seismograms,  and  real  earthquake  records. 

For  earthquake  seismograms,  we  visualize  a  choice  of  primaries 
made  by  the  criterion  discussed  above  as  implying  that  the  earth  pro¬ 
duces  the  minimum  number  of  primary  events.  Thus  the  function  r(t) 
in  equation  (1)  will  be  representative  of  an  earth  with  the  fewest 
number  of  major  discontinuities  as  opposed  to  one  with  many  minor 
boundaries.  This  implies  that  the  mot  ,  correct  set  of  primaries  rela¬ 
tive  to  the  other  inverse  filtered  seismograms,  will  concentrate  its 
energy  in  the  fewest  number  of  large  pulses.  Or  in  other  words,  it 
will  be  visually  the  simplest  in  structure.  For  this  reason  we  refer 


to  the  criterion  as  the  simplicity  criterion. 
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It  is  useful  for  illustrative  purposes  to  plot  C  as  a  function  of 
delay  time  T  when  u  is  held  constant.  The  value  of  T  for  which  this 
function  achieves  a  maximum  when  u  is  a  minimum  provides  the  estimate 

A.  i»>, 

T  for  the  correct  value  of  delav  time.  T  will  be  the  value  of  T,  used 

k 

A 

to  produce  the  estimated  primaries,  P.  In  the  next  section  we  discuss 
how  T  is  converted  into  focal  depth. 


Conversion  of  T  into  Focal  Depth:  The  seismogram  itself  cannot 
give  us  a  direct  estimate  of  the  focal  depth.  However,  if  the  velocity 
structure  in  the  area  of  the  source  is  known  and  the  value  of  T  is  ob¬ 
tained  from  the  seismogram,  then  the  focal  depth  can  be  computed  with 
travel-time  equations.  However,  we  cannot  expect  this  velocity  in¬ 
formation  to  be  available  in  the  area  of  each  seismic  disturbance. 

Since  this  is  the  case,  it  is  suggested  thaf  an  average  velocity 
structure  could  be  constructed  from  known  P-wave  velocities  as  a 
function  of  depth  .  (Birch  et  al  (1942)  and  Gutenberg  (1955,  1958) 
have  compiled  this  information).  This  average  velocity  profile  may 
be  mathematically  approximated  by  a  piecewise  linear  velocity  function 
(Steinhart,  1961)  of  the  form 

V.  ,  =  V.  f  b.(Z.  -  -  Z.)  (18) 

J+1  J  J  J+1  J 

where  V.  is  the  F-wave  velocity  at  depth  Z.,  b.  is  the  rate  of  change 
J  J  J 

tH 

of  velocity  with  depth,  and  the  subscript  "j"  refers  to  the  j  second 
order  velocity  discontinuity  in  the  downward  direction.  We  can  com¬ 
pute  the  time  for  a  P-wave  to  travel  from  Z .  to  Z 


using 


L  V-"  V 


(Steinhart,  1961).  Here  cp^  is  the  angle  of  incidence  of  the  P-wave 
path  of  travel  with  a  discontinuity  at  depth  Zy  The  horizontal  dis¬ 
tance  traversed  by  the  P-wave  going  from  Z.  to  Z .  ,  is  given  bv 

O  ]  O 


<.  =  r~*^ - : - -  (cos  tp.  -  cos  ip.  .). 

J  bj  sin  <P^  ^.1  j+i 


The  angle  of  incidence  of  P-waves  with  the  earth's  surface  for 
various  distances  has  been  theoretically  computed  by  Jeffreys  and 
Bullen  (1940).  Nuttli  and  Wliitmore  (1961)  made  an  observational  de¬ 
termination  of  P-wave  incidence  angles.  This  information,  combined 
with  observed  travel  times  for  a  zero  focal  depth  P-wave  (Jeffreys, 
1940)  and  the  equations  (19)  and  (20)  will  enable  us  to  construct 
travel-time  curves  for  P  and  pP,  for  different  focal  depths.  The  pro¬ 
cedure  for  doing  this  is  as  follows;  First  a  depth  of  focus  Z.  and  an 

J 

epicentral  distance  d  are  chosen.  For  this  distance  the  angle  of  in¬ 
cidence  of  a  P-wave  is  found.  Then  t.  and  x.  are  computed  for  the 

J  J 

depth  Zy  To  get  the  time  of  travel  of  the  ghost  (pP),  t^  is  added 
to  the  P-wave  zero  focal  depth  travel  time  (T^)  for  the  distance  d. 

The  epicentral  distance  for  pP  for  the  focal  depth  v/ili  be  d  ^  x.. 
Now  for  the  P-wave  generated  at  this  depth  the  reverse  is  true.  Tne 


time  t.  is  subtracted  from  T  and  its  epicentral  distance  will  be 

d  -  Xj.  In  a  similar  manner  the  remainder  of  the  travel-time  curve  can 

be  constructed  for  Z.  by  varying  d.  The  parameters  involved  are 

J 
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illustrated  in  Figure  11. 


Fig.  11.  Parameters  involved  in  constructing  travel-time  curves 
for  P  and  pP. 

This  completes  the  procedure  to  be  evaluated  in  this  thesis.  It  has 
three  distinct  advantages: 

1.  No  estimate  of  earth  filtering  is  required. 

2.  A  direct  estimation  of  T  can  be  obtained. 

3.  Subjectivity  by  the  user  is  minimized. 

A  disadvantage  of  the  method  is  that  the  length  of  the  inverse  operator 
is  directly  proportional  to  the  estimated  focal  depth.  For  large  focal 
depths  this  requires  operating  on  a  longer  portion  of  the  input  signal. 
This  could  result  in  storage  problems  in  the  computer  and  also  increase 
the  computation  time  quite  rapidly. 

Scope  and  Limitation  of  Study 

This  study  systematically  evaluates  a  method  designed  to  estimate 
two  parameters  contained  in  a  banddimited  signal  that  obeys  the 


doublet  hypothesis. 

Synthetic  data  is  used  for  all  experimental  evaluation.  Earth¬ 
quakes  were  analyzed  to  test  both  the  procedure  and  the  assumptions 
on  "live"  (actual)  data. 

Only  the  most  simple  ghosting  filter  is  assumed.  However,  the 
technique  presented  can  easily  be  modified  to  accoimiodate  more  complex 
ghosting  situations,  with  a  corresponding  increase  in  computation. 


II 


PROCEDURE  OF  THE  INVESTIGATION 

Design  of  the  Optimum  Inverse  Filters 

Inverse  filter  arrays  were  generated  using  the  computer  for 
three  noise  to  signal  ratios  (R)  (low,  medium  and  high  noise  level) 
and  twelve  ghost  amplitudes  (R^)  ranging  from  0.1  to  1.2.  Each 
operator  contains  twenty  points.  This  number  of  points  was  chosen  as 
a  compromise  between  expense  of  running  the  computer  and  increased 
effectiveness  of  a  long  inverse  operator.  An  inverse  operate-  designed 
for  a  specific  R  and  combination  may  be  modified  for  any  delay  time 
T,  simply  by  making  the  spacing  between  filter  points  equal  to  that  of 
the  desired  delay  time.  This  is  true  because  the  inverse  filter  is  a 
sampled  function  (sampled  at  0.1  second  intervals)  which  is  zero 
everywhere  except  at  multiples  of  T  (Foster  et  al,  1964).  Appendix  II 
contains  diagrams  of  the  filter  arrays  used  in  this  thesis. 

It  is  interesting  to  observe  that  the  inverse  operator  resigned 
for  small  ghost  amplitudes  consists  essentially  of  one  point  and  hence 
the  convolution  operation  will  simply  return  the  original  signal.  The 
envelope  of  these  inverse  operators  is  another  interesting  feature. 

When  the  noise  to  signal  power  ratio  is  large  the  operator  effectively 
reduces  to  two  points,  one  positive  and  the  other  negative.  Such  an 
opereto’"  will  do  very  little  to  estimate  the  primaries. 

Both  these  features  are  to  be  expected,  since  for  low  amplitude 
ghosts  very  little  filtering  is  required  to  suppress  them.  Also,  when 


signals  are  recorded  with  a  high  background  noise  level  very  liiile  can 
be  done  to  recover  the  primary  signal. 

Synthetic  Signals 

For  the  purposes  of  evaluating  the  effectiveness  of  the  pro¬ 
cedure,  synthesized  synthetic  data  containing  three  primary-ghost  pairs 
were  generated  in  the  computer.  Important  signal  parameters  were  built 
into  the  data.  These  included  variable  delay  times,  ghost  amplitudes 
and  random  noise.  The  most  important  single  parameter  is  the  delay 
time,  since  the  effectiveness  of  the  procedure  is  based  on  its  ability 
to  estimate  T  for  small  primary  ghost  separations  in  the  presence  of 
noise  and  limited  bandwidth  of  the  signal. 

A  reasonable  crustal  model  was  assumed  that  would  generate  at 
least  three  primary-ghost  pairs  in  a  short  interval  of  time.  Figure 
12a  is  a  diagram  of  the  crustal  layering.  It  represents  an  average  of 
several  profiles  given  by  Steinhart  (1961).  Three  r(t)  functions  were 
derived  for  focal  depths  of  2.5,  5.0  and  10.0  kilometers.  The  ampli¬ 
tudes  cf  the  primaries,  were  determined  from  amplitude  versus  angle  of 
incidence  curves  based  on  Zoeppritz  equations  (Steinhart,  1961). 

Figure  12b  illustrates  the  three  r(t)  functions  after  convolution  with 
h(t).  A  "modified"  Benioff  Geneva-Type  seismograph  system  impulse  re¬ 
sponse  b(t)  was  used  to  synthesize  seismograms  by  convolution  with  the 
h(t)  filtered  r(t)  function.  By  "modified"  we  mean  that  the  impulse 
response  used  in  this  thesis  has  been  expanded  and  is  given  by  b(t/2) 
rather  than  b(t).  This  results  in  shifting  the  spectrum  of  the 
synthetics  towards  lower  irequencies;  an  attempt  to  simulate 
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Fig*  12a.  Crustal  layering  used  to  generate  synthetic  spikograas 
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Fig.  12b,  Synthetic  spikograas  for  three  focal  depths 
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attenuation  in  the  earth.  The  impulse  I'esponse  curve  and  its  trequency 
spectrum  are  illustrated  in  Figure  13.  Random  noise  was  generated  in 
the  computer  with  noise  to  signal  power  ratio  equal  to  0.01,  0.0625  and 
0.25,  These  noises  were  then  added  to  the  filtered  r(t)  function  to 
produce  the  synthetic  seismograms.  The  synthetic  signals  have  the 
form 

s(t)  =  b(t)*h( t)*r( t)  +  n(t). 

By  combining  three  groups  of  variables,  T,  and  R  in  all 
possible  independent  combinations,  twenty-seven  different  signals  were 
generated.  Table  1  lists  the  values  for  these  parameters. 


TABLE  1.  Values  for  T,  R^  and  R  used  to  construct  synthetic 
seismograms . 


Focal  depth  in  km 

T  (sec) 

Ro 

R 

2.5 

0.3 

0.4 

0.0100 

5.0 

0.9 

0.7 

0.0625 

10.0 

1.6 

1.0 

0.2500 

Figure  14  is  a  block  diagram  of  the  evaluation  procedure  for 
S3mthetic  signals.  The  first  step  is  the  selection  of  three  inverse 
filters  designed  for  the  R  built  into  the  signal  s(t).  Each  one  is 
designed  for  a  different  ghost  amplitude.  A  range  of  delay  times  is 
chosen.  Only  one  ghost  amplitude  and  delay  time  combination  is  correct 
for  s(t).  The  signal  is  then  convolved  with  each  filter  designed  for 


procedure 


every  delay  time  and  ghost  amplitude  combination  possible.  This 
amounts  to  thirty  separate  convolutions.  The  simplicity  of  each  output 
signal  is  measured  for  ten  different  u-levcls  (eg.  10  to  100%  of  ma.xi- 
nium  signal  amplitude  in  10%  increments)  using  the  simplicity  criterion. 

A  A 

The  set  of  primarie.s  P(t;R^,T)  that  most  nearly  satisfies  the 
criterion  is  plotted  together  with  s(t)  and  the  appropriate  inverse 

A  A 

filter.  R  and  T  are  tlie  estimated  para.meters.  If  two  or  more  output 
o 

signals  appear  to  satisfy  the  criterion  equally  well,  then  they  are 
all  plotted,  so  that  they  may  be  subjectively  examined.  This  whole 
operation  is  performed  on  a  computer.  Appendix  III  describes  the  two 
main  program.s  and  several  subroutines  writcen  to  evaluate  the  delay 
time  resolution  for  the  procedure. 

Efi.  thquake  seismograms  were  used  to  further  evaluate  the  h  pth 
of  focus  estimation  procedure.  Theoretically  determined  delay  times 
were  converted  into  depths  in  order  that  they  could  be  compared  with 
published  focU  depths  for  the  same  earthruakes.  Since  very  little 
published  Information  is  available  for  focal  depths  less  than  33  kilo¬ 
meters,  it  was  necessary  to  construct  travel-time  curves  according  to 
the  procedure  discussed  on  pages  24,25,26.  Tic  lack  of  published 
information  for  shallow  focal  depths  further  reflects  the  failure  of 
visual  recognition  techniques  to  obtain  delay  times  for  reasons  already 
discussed,  such  as  overlap  of  pulses  and  high  background  noise  level. 

Construction  of  the  Travel-Time  Curves 

A  velocity  stinjcture  for  the  crust  and  mantle  was  construrted 


,;rom  available  data.  This  structure  combined  with  travel -times  for  a 


zero  focal  depcii  P-pulse  obtained  from  the  Jet  treys-Bul  len  SeisiTiologi - 
cal  Tables  (1940)  was  used  to  calculate  travel -times  for  P  and  pP,  Ihc 
calculations  v/ere  perfumed  by  a  computer  for  dept^is  of  focus  ranging 
from  5  to  100  kilometers,  and  epicentral  distances  from  500  to  10,000 
kilometers.  One  of  the  more  useful  plots  of  the  travel -times  may  be 
seen  in  Figure  15.  Here  depth  of  focus  is  plotted  as  a  function  of 
the  travel-time  difference  of  P  and  pP  for  constant  epicentral 
distances.  It  is  readily  apparent  that  the  time  difference  is  almost 
independent  of  distance  for  depths  less  than  20  kilometers.  Further, 
it  is  apparent  that  the  rate  of  change  of  delay  time  with  epicentral 
distance  decreases  as  a  function  of  epicentral  distance  for  constant 
focal  dep^h.  The  results  of  this  theoretical  model  agree  quite  well 
with  available  data  from  earthquakes  with  depths  of  focus  greater  than 


50  kilometers. 


differenc*  darivcd 


EXPERIMENTAL  RESULTS 


Single  Doublet  Function 

As  an  initial  test  the  procedure  was  evaluated  using  a  single 
doublet  function  (Figure  16  bottom  diagram).  The  negative  impulse 
was  four  tenths  the  amplitude  i  f  the  positive  spike.  The  delay  time 
was  half  a  second  corresponding  to  four  sample  points  in  the  function. 
Random  noise  at  e.  level  of  0,1  relative  to  the  doublet  was  added. 

For  the  analysis  three  inverse  filters  were  designed  for  ghost 
amplitudes  of  0,2,  0,4  and  0.6.  The  noise  to  signal  power  ratio  used 
in  the  derivation  of  the  inverse  filters  was  0.01.  Delay  times  varied 
from  0.]  tc  1.1  seconds  in  increments  of  0.1  seconds.  Each  filter  was 
modified  for  every  delay  time  in  the  range  and  convolved  with  the 
doublet.  The  deconvolved  signal  was  then  measured  ^or  structural 
simplicity  for  ten  different  u-levels. 

The  results  of  this  analysis  are  presented  in  Figures  16a  through 
16d.  Figures  16a  and  16b  are  plots  of  signal  simplicity  (energy  con¬ 
centration)  as  a  function  of  delay  time  for  specific  u-levels.  For 
u  =  0.1  there  is  2  definite  high  energy  concentration  corresponding  to 
the  filter  designed  for  the  correct  delay  time.  The  largest  peak 
corresponds  to  the  filter  that  was  also  designed  for  the  correct  ghost 
amplitude. 

For  the  u-level  0.5  the  magnitude  of  the  peaks  are  considerably 
reduced.  However,  they  do  retain  the  same  spacial  orientation  ‘ith 
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respect  to  T.  Both  sets  of  curves  indicate  that  our  theoretical  con¬ 
siderations  are  correct. 

To  further  demonstrate  the  ability  of  the  procedure  for  selecting 
the  correct  delay  time  and  ghost  amplitude,  two  more  sets  of  energy 
concentration  curves  were  constructed.  Figure  16c  is  a  plot  of  energy 
concentration  for  the  correct  delay  time  versus  u-level.  For  small 
u-levels  the  curves  are  separated  and  indicate  that  the  procedure  pro¬ 
duces  a  maximum  for  the  correct  ghost  amplitude  as  expected.  However, 
for  higher  u-levels  the  energy  curves  tend  to  become  superimposed. 

This  is  not  the  case  when  we  consider  energy  curves  for  filters  designed 
for  the  correct  ghost  amplitude  but  wrong  delay  time.  Figure  16d 
clearly  illustrates  this  fact. 

Based  on  this  analysis  it  is  concluded  that,  the  procedure  is 

_  f 

able  to  select  from  an  array  of  inverse  filters  thdtytilter  which  is 
designed  for  the  ghost  amplitude  and  delay  time  built  into  the  doublet 
function. 

\ 

BandTimited  Single  Doublet  Functions 

Synthesized  doublet  functions  were  constructed  by  convolving  the 

! 

same  doublet  with  the  modified  Benioff  impulse  response  illustrated  in 
Figure  13.  Three  delay  times  were  chosen  so  that  the  following 
conditions  would  occur.  First,  for  a  delay  time  of  2,1  seconds  the 
primary  pulse  would  be  completely  separated  from  its  ghost  reflection 
Second,  for  a  pulse  separation  of  0,9  seconds  they  would  partially 
overlap.  And  third,  for  a  pulse  separation  of  0.5  seconds  overlap 
would  occur  to  such  a  degree  that  individual  pulses  would  be  visually 


difficult  to  identify. 


Each  signal  waa  tested  in  the  same  manner  that  the  impulse 
doublet  was  analyzed.  Figure  17  shows  the  results.  For  the  signal 
with  a  primary-ghost  delay  time  corresponding  to  0.5  seconds,  sub¬ 
sidiary  maxima  appear  in  the  energy  concentration  versus  delay  time 
curve.  However,  the  maximum  corresponding  to  the  correct  delay  time 
is  the  most  pronounced.  It  should  be  noticed  by  comparison  with 
Figures  16a  and  16b  that  band-limiting  does  decrease  the  sensitivity 
(measured  in  terms  of  peak  amplitudes)  of  the  procedure. 

Using  the  location  of  maxima  as  a  guide,  a  filter  was  chosen 
which,  convolved  with  the  input  signal,  generated  the  best  estimate 
of  the  primaries.  The  results  of  this  operation  for  the  three 
synthesized  doublets  is  shown  in  the  lower  half  of  Figure  17.  In  all 
cases  the  original  pulse  shape  has  returned  (compare  with  impulse  re¬ 
sponse  Figure  13),  and  the  ghost  has  been  effectively  suppressed. 

Thus  for  these  three  cases  the  procedure  has  been  successful  in  es¬ 
tablishing  which  combination  of  and  T^^  was  correct. 

Bandiimlted  Triple  Doublet  Functions 

Before  testing  the  criterion  on  earthquake  seismograms  a  sequence 
of  synthetic  signals  was  generated  that  contained  three  primary-ghost 
combinations.  The  construction  of  these  signals  is  discussed  in  Part 
11  (Synthetic  signals).  Three  primary  pulses  were  chosen  to 
simulate  actual  earthquakes  with  shallow  focal  depths  at  small  epi- 
central  distances.  Small  epicentral  distances  were  considered  to  Be 
a  greater  challenge  to  the  procedure  for  the  same  number  of  primary 
pulses,  since  proximity  of  primaries  is  more  likely  than  for 


Simplicity  criterion  uaed  to  determine  the  correct  delay  time  for  three 

band-limited  doublet  functione  containing  random  noiae. 


teleselsmic  (large)  epicentral  distances.  For  teleseismic  distances, 
the  procedure  would  yield  essentially  the  same  results  as  in  the  case 
of  the  single  bandiimited  doublet,  since  the  synthetic  seismograms 
would  essentially  consist  of  a  set  of  separated  primary-ghost  pairs. 

Synthetic  Signals  for  a  2.5  km  Focal  Depth  Nine  synthetic 
signals  for  a  delay  time  of  0.3  seconds  were  analyzed.  Figure  illc- 
strates  these,  together  with  the  deghosted  signals  selected  by  ttic 
simplicity  criterion  (Figure  19).  Several  features  in  these  figures 
are  worth  emphasizing. 

With  reference  t.  Figure  18.,  synthetic  signals  for  large  ghost 

amplitudes  (R  =  1.0)  are  quite  oscillatory.  Individual  pulse  recog- 
o 

nition  is  virtually  impossible.  In  these  cases  the  inverse  operator 
designed  for  the  correct  ghost  filter  failed  to  improve  the  visual 
character  of  the  synthetic  signal. 

The  deghosted  signal  P(t)  is  in  general  broader-band  than  the 
Input  signal  S(t).  Signals  for  =  1.0,  R  =  0,01  and  R^  =  0.7,  R  = 
0.25  particularly  emphasize  this. 

For  R  =  0.7,  R  =  0.25  and  R  =  1.0,  R  =  0.625  and  R  =  0.25, 
o  o 

changing  the  desigr.  of  the  inverse  operator  did  not  substantially 
effect  the^energy  concentration  versus  delay  time  curves.  This  indi¬ 
cates  that  the  procedure  seems  to  lose  sensitivity  as  the  ghost  ampli¬ 
tude  and  noise  to  signal  power  ratio  increase. 

Except  for  R^  =  1.0,  R  =  0.01,  the  procedure  selected  the  correct 
delay  time  built  into  these  nine  synthetic  signals.  For  this  exception 
a  delay  time  of  0.2  seconds  was  chosen.  This  is  only  0.1  second  less 
than  the  correct  delay  time  (0.3). 
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Figj  l8.  Synthetic  signals  SCt)  for  a  2.5  km  focal  depth 

(0.3  sec  delay  time)  with  their  deghosted  signals  P(t) 
selected  by  the  simplicity  criterion. 
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Synthetic  Signals  for  a  S.O  km  Focal  Depth,  Nine  synthetic  sig¬ 
nals  for  a  delay  time  corresporic ing  to  0,9  seconi’s  were  analyzed. 

Figure  20  illustrates  these  together  with  the  deghosted  signals  se¬ 
lected  by  the  simplicity  criterion. 

Curves  of  energy  concentration  versus  delay  time  for  each  signal 
are  illustrated  in  Figure  21.  For  several  cases  these  curves  have  two 
ma.xima,  one  at  0.3  seconds  and  the  other  at  0.9  seconds.  The  maximum 
at  0.3  seconds,  although  pronounced,  is  for  an  incorrect  inverse  filter 
which  when  convolved  with  the  input  signal  fails  to  simplify  it  (in 
the  sense  defined  in  Primary-estimation  Method).  The  cases  where  this 

A  A 

is  most  obvious  are  R  =  0.6,  R  =  0.0625  and  R  -0.5,  R-  0,25.  For 

0  o 

these  cases  the  simplicity  criterion  was  unable  to  differentiate  the 
correct  from  the  incorrect  inverse  filtered  record.  The  reason  for 
this  seems  to  be  that  energy  concentrat ion  increases  with  increasing 
bandwidth.  Thus  if  an  incorrect  inverse  filter  overcompensates  the 
spectrum  of  the  input  signal,  as  in  the  discussion  on  page  19  ,  then 
the  simpllcicy  criterion  may  choose  the  incorrect  version  by  virtue  of 
its  broader  bandwidth  over  the  correct  version.  This  is  a  limitation 
of  the  procedure. 

The  procedure  definitely  helped  to  select  the  correct  delay  time 

in  5  out  of  the  9  synthetic  signals  analyzed.  For  R^  0.4  and  R  - 

0.25  the  procedure  selected  1.0  and  0,8  seconds  for  1.  Either  of 

these  values  differ  only  by  0. 1  second  from  the  correct  value  of  I 

(0,9).  For  the  three  noise  levels  with  R  =1.0  the  procedure  was 

o 

unable  to  estimate  any  delay  time. 
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Fig,  20,  Synthetic  signals  S(t)  for  a  5.0  km  focal  depth 

(0.9  sec  delay  time)  with  their  deghosted  signals  P(t) 
selected  hy  the  simplicity  criterion. 


concentration 
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Synthetic  Signals  for  fi  IQ.Q  km  Focal  Depth:  Nine  synthetic 
signals  for  a  delay  time  of  1.6  seconds  were  analyzed.  Figure  22  illu¬ 
strates  these  together  witn  the  deghosted  signals  selected  by  the 
simplicity  criterion. 

The  set  of  energy  concentration  curves  for  this  delay  time  are 

uniform  in  visual  character  (Figure  23).  Only  one  maximum  appears  in 

each  sot  of  curves  at  a  delay  time  of  1.6  seconds.  These  maxima  give 

correct  values  for  both  T  and  R  in  all  cases,  except  for  R  =:  0.4, 

o  o 

R  =  0.0625;  R  0.4,  R  =  0. 25;  and  R  =  1.0,  R  _  0.25.  For  these 
o  o 

three  cases  the  procedure  selected  the  correct  delay  time  built  into 
the  synthetic  signal,  but  not  the  correct  R^. 

The  second  ghost  (third  pulse)  was  not  suppressed  when  the 
synthetic  signal  Si(t)  contained  ghosts  with  amplitudes  greater  than 
0.7,  Even  though  this  may  have  been  the  case,  the  procedure  was  still 
able  to  select  the  correct  delay  time. 

Analysis  of  Earthduake  Seismograms 

Four  earthquake  seismograms  were  chosen  to  test  the  estimation 
procedure  on  the  basis  of  visual  structure  and  the  depth  of  focus  as 
published  by  the  USC6tGS.  Each  earthquake  was  assumed  to  obey  our 
mathematical  mode]  .  A  delay  time  and  a  ghost  air  litude  was  found  for 
each  one  using  the  procedure.  The  delay  time  was  converted  into  focal 
depth  using  travel -time  curves  constructed  from  a  piecewise  continuous 
velocity  function. 

The  deepest  earthquake  analyzed  occurred  in  the  Kurile  Island 
region  at  a  depth  of  approximately  60  kilometers.  Its  visual  structure 


Fig.  22.  Synthetic  signals  S(t)  fo:  -  i;.0  km  focal  depth  • 

(1.6  sec  delay  time)  with  their  deg^osted  signals  P(t) 
selected  by  the  simplicity  criterion. 


oynthatic  signals  generated  for 
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is  relatively  unconipl icated  (Figure  24  top).  Apart  from  apparent  noise 
only  two  large  pulses  are  in  the  interval  analyzed,  A  series  of  in¬ 
verse  operators  were  applied  to  the  earthquake  assviming  a  noise  to 
signal  power  ratio  of  0.25.  This  assumption  appeared  reasonable  from 
visual  inspection.  Application  of  the  procedure  strongly  indicated 
that  the  delay  time  is  12.7  seconds  (Figure  25),  which  according  to  our 
travel-time  curves  corresponds  to  45  kilometers,  A  ghost  amplitude  of 
1.0  was  estimated  by  the  procedure.  After  applying  the  operator  de¬ 
signed  for  these  values  a  seismogram  was  obtained  which  is  very  similar 
to  the  input  seismogram,  except  that  the  second  pulse  is  completely 
suppressed.  It  should  be  noticed  that  the  original  pulse  shape  of  the 
primary  pulse  has  not  been  altered  after  deconvolution  (Figure  24). 

The  second  earthquake  analyzed  has  a  published  focal  depth  of  44 
kilometers.  It  occurred  near  Co.sta  Rica  almost  3500  kilometers  from 
State  College,  Pennsylvania.  An  analysis  of  tnis  signal  produced  two 
prominent  peaks  in  the  energy  concentration  versus  delay  tim.e  curve 
for  a  ghost  amplitude  of  0.8  (Figure  25).  Two  deconvolved  seismograms 
were  obtained  using  both  the  operators  indicated  to  be  most  correct 
by  the  procedure.  The  operator  designed  for  a  delay  time  corresponding 
to  6.9  seconds  has  slight  y  suppressed  the  pulses  following  tbe  primary, 
but  has  also  boosted  the  apparent  noise  level.  The  second  operator 
designed  for  a  delay  time  of  8.7  seconds  has  successfully  suppressed 
two  pulses,  but  boosted  the  ^hird  pulse.  This  pulse  has  a  shape  almost 
identical  in  character  to  the  first  primary  pulse.  It  appears  too 
early  in  time  to  be  associated  with  any  known  arrival  for  this  e.arth- 
quake.  It  has  been  suggested  that  if  the  8.7  second  operator  is 
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Fig,  24,  Input  earthquake  seiemograms  and  their  deghosted 
counterpart  selected  by  the  simplicity  criterion. 
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correct,  then  this  pulse  is  the  first  primary  pulse  from  a  second 
earthquake  at  the  same  source.  Any  conclusions  in  the  matter  dra™ 
from  a  single  earthquake  seismogram  represent  only  a  judgement.  It 
should  be  emphasized  that  the  criterion  cannot  distinguish  between  the 
two  operators  from  the  information  given.  However,  the  fact  that  the 
apparent  noise  is  boosted  when  the  shorter  operator  is  used  makes  one 
sujpicious  of  this  delay  time.  The  longer  delay  time  corresponds  to  a 
depth  of  focus  of  31.5  kilometers. 

Analysis  of  the  third  earthquake,  whose  source  was  in  Alaska  more 
than  5000  kilometers  from  State  College,  indicated  that  its  depth  of 
focus  was  10  kilometers.  Its  published  focal  depth  is  25  kilometers. 
This  earthquake  is  relatively  simple  in  structure  with  only  two  major 
pulses.  The  second  pulse,  if  it  is  a  ghost,  is  more  ocillatory  than 
the  primary.  The  delay  time  for  the  maximum  in  the  energy  concen¬ 
tration  curve  was  3.3  seconds  (Figure  24).  Other  maxima  appear  in  the 
curve  but  the  criterion  dictates  that  the  correct  operator  yields  the 
greatest  concentration  of  energy  above  the  minimum  value  of  u  for  the 
most  likely  delay  time.  This  is  the  operator  that  should  be  chosen, 

The  inverse  filter  corresponding  to  the  3.3  second  delay  time  succ;  so- 
fully  suppressed  the  pulse  following  the  primary.  The  rest  of  the 
Signal  remains  relatively  unchanged  (Figure  24), 


The  last  earthquake  analyzed  occurred  at  the  Hexic 


u-uurfuemaia 


border,  and  was  computed  by  the  USC4G3  to  have  a  12  kilometer  focal 

depth.  The  procedure  estimated  T  =.  4.9  sec.  and  R  =  1.0  1  =  49 

o 

sec.  corresponds  to  a  focal  depth  of  15  kilcmeter.s.  This  shaJlcv 
earthquake  seismogram,  bcs  a  v.i3uaily  complicated  structure  (Figure  2^. 
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last  earthquake  pair).  Indi.’idual  pulses  cannot  be  easily  identified. 
Hov,ever,  the  procedure  ”cry  strongly  suggests  what  and  T  should  be. 
There  is  no  ambiguity  at  all  in  this  case  (Figure  25).  The  decon- 
volved  seismogram  (Figure  24)  is  substantial ly  simpler  in  structure 
compared  to  the  original  seismogram.  Most  of  the  pulses  following  the 
primary  have  been  suppressed  and  no  additional  pulses  have  been 
created . 


SUMMARY  ANL  CONCLUSIONS 


This  thesis  is  concerned  with  the  problem  of  extracting  infor¬ 
mation  from  a  band-limited  signal  in  the  presence  of  corrupting  noise. 

In  particular  the  signal  is  a  seismic  signal  generated  at  or  near  the 
earth's  surface.  The  information  desired  is  a  correct  estimate  of 
focal  depth  of  the  earthquake's  source. 

An  accurate  estimate  of  the  source  depth  is  a  criterion  that 
would  help  to  differentiate  between  natural  earthquakes  and  clandestine 
nuclear  explosions, 

Currently  the  USC&GS  routinely  publishes  the  focal  depth  of  most 
earthquakes  recorded  on  a  world-wide  basis.  This  depth,  especially  if 
the  earthquake  is  shallow,  oiay  be  in  rror  by  as  much  as  2b  kilometers 
(Gunst  and  Engdahl,  1962). 

Several  investigators  have  suggested  techniques  that  will  give 
focal  depth  information.  These  techniques  fall  into  three  main 
categories:  visual  recognition  procedures,  spectral  analyses  and 
linear  filter  operations. 

This  thesis  evaluates  a  technique  for  extracting  delay  time  in¬ 
formation  from  a  signal  that  obeys  a  specific  mathematical  model. 

This  model  is  assumed  to  be  a  reasonable  approximation  of  the  actual 
data  generation  process  in  the  earth. 

The  method  is  based  on  the  assumption  that  every  seismic  signal 
has  passed  through  a  ghost  filter.  This  filter  introduces  a  pulse 
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(gb.ost)  after  each  primary  pulse  in  the  signal.  It  is  delayed  by  a 
time  cliat  is  related  to  the  focal  depth  of  the  earthquake.  Further, 
it  i?  assumed  that  its  polarity  is  reversed,  but  it  retains  the  primary 
pulse  .shape.  The  ratio  of  the  amplitude  of  the  primary  to  its  ghost  is 
permitted  to  vary. 

The  estimation  procedure  uses  optimum  inverse  filters  of  the 
doublet  function  (ghost  filter)  together  with  a  criterion  that  measures 
the  visual  simplicity  of  a  seismic  signal  convolved  with  an  inverse 
filter.  The  inverse  filters  are  designed  to  extract  primary  energy  in 
the  presence  of  ghost  energy  and  random  noise  on  the  basis  of  a  i  .-asc- 
mean-square  error  criterion.  Filter  design  is  dependent  on  the  delay 
time  betv/een  primaries  and  ghosts,  their  amplitude  ratio  and  noise  to 
signal  power  ratio.  The  criterion  vjas  devised  on  the  basis  that  a 
maximum  in  the  energy  concentration  of  the  inverse  filtered  seismograms 
above  a  minimum  level  indicates  tvhich  filter  was  most  correctly  de¬ 
signed.  This  Is  visualized  as  an  expression  of  the  assumption  that  the 
primary  seismogram  is  generated  by  a  few  large  discontinuities,  rather 
than  by  many  minor  boundaries. 

The  technique  was  applied  to  bandTimited  synthetic  signals  that 
contained  several  primary-ghost  pairs  in  the  presence  of  random  noise. 
Each  signal  was  convolved  v/ith  an  array  of  inverse  operators,  one  of 
which  was  designed  for  the  delay  time,  amplitude  ratio  and  noise  level 
built  into  the  signal.  The  criterion  was  then  applied  to  determine 
which  operator  was  most  correct.  Of  27  synthetic  signals  which  v/ere 
analyzed,  the  procedure  successfully  selected  the  correct  delay  time 
in  22  cases.  Figure  26  contains  plots  of  the  results  of  the  analyses 
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Fif,  26.  Graj^  aQaaa  'iaisje;  the  T,  R  aelectlon  capability  of  the 
eatiaatiwi  procedure  for  the  aynthetic  signals. 


for  three  different  noise  to  signal  power  ratios. 

Four  actual  earthquake  seismograms  were  then  atialyacd  in  the  same 
manner.  Focal  depths  computed  from  the  selected  delay  times  appeared 
quite  reasonable  when  compared  with  depths  for  the  same  earthquakes 
published  by  the  USC&GS.  In  all  cases  the  inverse  operator  selected  by 
the  criterion  considerably  simplified  the  original  signal. 

On  the  basis  of  the  results  of  the  analyses  of  synthetic  and 
earthquake  signals  it  is  concluded  that  the  method  has  several  distinct 
advantages.  The  three  most  important  advantages  are  listed  below. 

1.  No  estimate  of  earth  filtering  is  required. 

2.  A  direct  estimate  of  the  delay  time  between  pri¬ 
maries  and  ghosts  is  obtained. 

3.  Subjectivity  by  the  user  is  minimized. 

The  following  five  statements  summarize  the  effectiveness  of  the  pro¬ 
posed  focal  depth  estimation  procedure. 

1.  Analysis  of  synthetic  signals  indicates  that 
the  sensitivity  of  the  criterion  decreases  with 
increasing  ghost  amplitude,  a  decrease  in  band¬ 
width,  and  large  noise  to  signal  power  ratio. 

2.  If  more  than  one  filter  is  predicted  to  be 
equally  likely,  then  a  visual  examination  of  the 
dc'^onvolved  signals  may  be  sufficient  to 
determine  which  is  correct. 

3.  The  inverse  filter  selected  by  the  criterion 
docs  simplify  the  visual  character  of  the  input 
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4.  The  criterion  in  most  cases  is  successful  in 
choosing  the  correct  inverse  operator  for 
band-limited  synthetic  signals  of  the  type  used. 

5.  Analysis  of  the  earthquake  seismograms  indi¬ 
cates  that  the  ghost  filter  assumption  may  be 
reasonable. 

The  method  may  be  refined  in  several  ways.  First  of  all,  the 
ghost  filter  may  be  made  more  complicated  to  account  for  additional 
pulses  generated  between  the  source  and  the  earth's  surface.  The 
inverse  operators  could  be  made  longer  and  hence  more  effective  in 
reducing  the  doublec  or  ghosting  function  into  a  single  impulse 
(Watson,  1964:  personal  communication).  The  criterion  could  be  modi¬ 
fied  to  remove  its  bias  towards  band-broadened  signals.  In  this  way 
some  ambiguity  might  be  removed  and  hence  less  subjectivity  would  be 
invo Ived. 

The  estimation  procedure  certainly  is  encouraf'ing  enough  to 
warrant  further  investigation.  To  establish  fully  the  effectiveness 
of  the  procedure  in  accurately  determining  earthquake  focal  depths, 
several  earthquakes  recorded  by  many  different  stations  should  be 
analyzed.  Consistency  in  the  results  obtained  would  also  be  a  strong 
indication  that  the  doublet  hypothesis  is  a  good  assumption  for  earth¬ 
quakes  . 


APPENDIX  A 


Proof  of  simplicity  criterion  for  doublet  function 
with  no  noise  and  perfect  inverse. 


and  n  Impulses 


jd/(t-Tp|  -uj-  .  fd^(t-Tj)|  -a,]2  |d„S(t-yl  -uJ2 

nht  f  d„^^(t-T2)^  -  ...  ^  d^^cf(t-T^)^] 

d,  ,  d„ ,  .  .  .  ,  d  >u  .  . 

12  n  1 

Now  for  u  -  0,  c(0;R  .,T,  ),  -  and  c(0;R  .,1.)  -  — ^  . 

oj  k  1  £st  oj  k  n  ni^ii: 

Hence  c(0;R  .,T,  ),  >  c(0;R  .,T,  )  . 

OJ  k  1  OJ  k  n 

It  is  obvious  that  c(0;R  .  ,T,  'i ,  >c(u .  ;R  .,T,  ),  for  u.>0. 

OJ  k  1  1'  oj  k' 1  i 

For  the  correct  inverse  operator  c(u.;R  .,T.  )  will  be  a  maximuin  for 

I  OJ-  k' 


the  smallest  u-level.  That  is  for  a  doublet  function  with  no  noise 
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the  correct  and  ideal  operator  yields  a  single  impulse  as  desired. 
Further,  for  this  filter  no  other  maxima  will  appear  in  the  energy 
concentration  versus  u-Uvel  curve.  However,  maxima  may  appear  in  this 
curve  for  incorrect  filters,  but  of  course,  not  at  the  minimum  u-level. 
This  argument  may  be  extended  to  more  than  one  doublet  with  similar 
conclusions . 

When  noise  is  introduced  into  the  system  the  minimum  level  for 
a  maximum  energy  concentration  is  raised.  The  arguments  of  the  proof, 
are  similar  to  those  outlined  for  the  noiseless  situation,  except  that 
an  actual  inverse  filter  is  used. 
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APPENDIX  B 


Optimum  inverse  operators  designed  for  the  doublet 


function. 
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o 


»0  ‘  0-3 


*o  =  O-** 


••••«•«  *4 


=  “’5 


«o  =  0-6 


r<  =  0.7 

o 


R  =  0.8 

o 


« 


»o  =  ^-9 


R  =1,0 
o 


R^  =  1.1 
o 


B  =  1.2 

o 


Fis.  28.  ^^..S  of  t«.l„  20-point  intor..  operator. 

aesignec*  for  the  doublet  function  for  R  =  0.01. 
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-c.cccc 

“C.CCCO 

-c.cccc 
—  C  *  C  u  c  c 
-c.cccc 

-C.CCCO 

-o.cccc 

—C.CCCC 

-C.CCOi 

-C.GCIC 

l.CCCO 

C.C99C 

C.CC98 

O.CCIC 

O.CCCI 

c.cccc 

c.cccc 

c.cccc 

C.CCCC 

c.cccc 

H  *0*1 
o 


-C.C311 

-C.C632 

-C.C973 

-C. 

-C.1757 
-C.2225 
l.CCCC 
C.835? 
C.697^ 
C.581C 
0.6826 
C.3994 
0. 3286 
C.26SI 
C.2i59 
C.  17C5 
C.  13  06 
C  .  CO  6  6 
C.C616 
C,C3C2 


fi  =0.7 

o 


Table  2. 


-C.CCCC 

-C.CCCC 

-C.CCCC 

-c.cccc 

-C.CCCC 

-c.cccc 

-c.cccc 

-c.cccc 

-C.CCCC 

-c.cccc 

-c.cccc 

-C.CCCC 

-C.CCCO 

-C.CCCl 

-C.CCCl 

-C.CCC3 

-C.CCC6 

-C.CC12 

-C.CC22 

-C.CC39 

l.CCCC 

l.CCCC 

r . 1979 

0.2965 

C.C391 

0.0879 

C.CC77 

C.C261 

C.CC15 

C.CC77 

C.CCC3 

C.0C23 

O.CCCl 

C.CCC7 

C.CCCC 

C.CCC2 

C.CCCC 

C.CCCl 

c.cccc 

C.CCCC 

R  =0.2 
o 

R^=0.3 

O 

-C.CC29 

-0.CC93 

-0.CC62 

-C.C192 

-C.01C5 

-0.C3C5 

-C.C166 

-0.0639 

-C.C266 

-C.C6C6 

-C.C366 

*  C.C812 

l.CCCC 

l.CCCC 

C.6&C6 

C.  7671 

C.6632 

0.5882 

C.3152 

C,65C9 

C.2165 

G. 3653 

C. 1659 

C.266C 

C.C992 

C.2CIA 

C.0676 

C,153C 

C.C656 

C.1153 

C.C317 

c.ce5e 

C.C2C6 

C.C626 

C.C131 

C.C633 

c.cc7e 

C,C2?3 

C.CC  36 

C.Cl 32 

R  =0.8  R  =0.9 
o  o 


-C.CCCC 

-C.CCCl 

-C.CCCC 

-0.0CC3 

-C.CCCl 

-C.CCC6 

-C.CCC2 

-C.CC13 

-0.CCC6 

-C.CC26 

-C.CC  1C 

-0.CC53 

-C.CC25 

-c.cice 

-0.CC65 

l.CCCC 

l.CCCC 

0.6918 

C. 3966 

C.2619 

C. 1557 

C. 1 19C 

C.C616 

0.0585 

0.0262 

c.c2ee 

C.CC96 

C.0162 

C.CC38 

C.CC7C 

C.CC  15 

C.CC36 

G.CCC6 

C.CC  17 

G.CCC2 

c.ccce 

C.CCCl 

C.CCC6 

C.CCCC 

C.CCC2 

R  =0,4 
o 

R  »0.5 
o 

-C.Q736 

-0.C317 

-C.1632 

-0. C666 

-C.226C 

-C.C989 

-C.3C83 

-C. 1362 

-0.3968 

-C. 1777 

-C.6932 

-0.2263 

-C.5995 

-C.2  776 

-0.7177 

-0.3387 

-C,85C6 

-G. 6C9e 

-l.CCCC 

-C.6929 

l.CCCC 

-C.59C3 

C.85Q6 

-C.7C5C 

0.7177 

-0.86C2 

C.5995 

-l.CCCC 

G.6932 

0.2667 

0.3968 

C.  1955 

C.3C83 

C. 1699 

C.226C 

0. 1C83 

C. 1602 

C.C7CB 

0.0736 

C.C269 

R  =1.0 
o 

R  =1.1 
o 

Values  for  twelve  20-point  inverse  operators 
designed  for  H  =  0.01, 


-o.cccg 

-0.CC2C 
-C.CC37 
-C.CC65 
-0.0112 
-0.0i9C 
l.CCCC 
0.5877 
C. 3653 
0.2C29 
U. 1193 
C.07C1 
C.0612 
C.C262 
C.C162 
C.CC83 
Q.0C6e 
C.CC2fi 
.GC15 
C.  CCC7 

R  »tO»6 
o 


-C.C166 
-0.C299 
-C.0669 
-C.C663 
“C.C896 
■  C.1172 
-0.1513 
-C.1936 
-C.2662 
-C.312C 
-C.3966 
-C.6986 
-C.62e9 
-0.7932 
-l.CCCC 
C.IC55 
C.C/83 
G.C553 
C.C353 
C.Ol 72 


R  =1.2 
o 


X 


-0.0000 

-0.0000 

“0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-O.OOOl 

-0.0000 

-O.OOOl 

-0.0005 

-O.OOOl 

-0.0005 

-0.0018 

-0.0006 

-0.0025 

-0.0063 

-0.0060 

-0.0131 

-0.0227 

l.OOOO 

l.OOOO 

l.OOOO 

0.0940 

0.  1873 

0.2792 

0.0086 

0.0351 

0.0780 

0.0008 

0.0066 

0.0218 

0.0001 

0.0012 

0.0061 

0.0000 

0.0002 

0.0017 

0.0000 

0.0000 

0.0005 

0.0000 

0.0000 

O.OOOl 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

0.0000 

R  »0.1 
o 

R  =0.2 
o 

R  =0.3 
o 

-0.0031 

-0.0060 

-0.0151 

-0.0070 

-0.0130 

-0.0323 

*0.0127 

-0.0223 

-0.0540 

-0.0218 

-0.0357 

-0.0830 

-0.0366 

-0.0554 

-0. 1234 

-0.0609 

-0.0849 

-0.1806 

-0.1012 

-0.1296 

-0.2626 

-0.1680 

-0.1973 

-0.3805 

1.0000 

-0.3002 

-0.5505 

0.6029 

1.0000 

1.0000 

0.3634 

0.6576 

0.6920 

0.2191 

0.4324 

0.4787 

0,1321 

0.2842 

0.3309 

0.0796 

0.1866 

0.2284 

0.0479 

0.1223 

0.1571 

0.0288 

0.0797 

0-1073 

0.0171 

0.0514 

0.0722 

O.OlOO 

0-0322 

0-0470 

0.0055 

0.0187 

0.0281 

0.0024 

0.0086 

0.0132 

R  *0.7 
o 

50 

o 

II 

o 

• 

00 

R  *0.9 
o 

-0.0000 

-0.0002 

-0.0009 

-0.0000 

-0.0005 

-0.0021 

-0.0000 

-O.OOll 

-0.0041 

-0.0001 

-0.0025 

-0.0079 

-0.0003 

-0.0056 

-0.0148 

-0.0007 

-0.0122 

-0.0279 

-0.0018 

-0.0270 

-0.0523 

-0.0050 

-0.0594 

-0.0980 

-0.0136 

1.0000 

l.OOOO 

-0.0368 

0.4541 

0.5333 

1.0000 

0.2062 

0.2844 

0.3687 

0.0936 

0-1517 

0.1359 

0.0425 

0.0809 

0.0501 

0.0193 

0.0431 

0.0185 

0.0088 

0.0230 

0.0068 

0.0040 

0.0123 

0.0025 

0.0018 

0.0065 

0.0009 

0.0008 

0.0034 

0.0003 

0.0004 

0.0017 

O.OOOl 

O.OOOl 

0.0007 

R  *0.4 

0 

R  *0.5 
o 

R  *0.6 
o 

-0.0213 

-0.0135 

-0.0103 

-0.0453 

-0.0288 

-0.0222 

-0.0750 

-0.0479 

-0.0377 

-0.1140 

-0.0735 

-0.0592 

-0.1673 

-0.1091 

-0.0902 

-0.2415 

-0.1593 

-0.1356 

-0.3459 

-0.2310 

-0.2029 

-0.4935 

-0.3337 

-0.3026 

-0.7028 

-0.4815 

-0.4509 

-l.OOOO 

-0.6940 

-0.6716 

l.OOOO 

-1.0000 

-l.OOOO 

0-7028 

0.5817 

0.3669 

0.4935 

0.4032 

0.2462 

0.3459 

0,2791 

0.1651 

0-2415 

0.1925 

0.1104 

0.1673 

0.1318 

0-0734 

0-1140 

0.0889 

0.0481 

0.0750 

0.0579 

0.0306 

0.0453 

0.0348 

0,0180 

0.0213 

0.0163 

0.0084 

R  =1.0 
o 

H  =1.1 
o 

R  *1,2 
0 

Table  3.  Values  for  tweiire  20-point  inverse  operators 
designed  for  8  =  0.0^25- 
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fi  =  0*1 
o 


R  =  0.2 

o 


8^  =  0.3 
o 


R  =  0.4 
o 


•••••••• 


•••••  M««*« 


**•••. 


%  »  0.5 


R  =  0.6 
o 


Ro  =  0.7 


R  =  0.8 
o 


*  i 


••  •*«  •••*, 


’•M* 


K  =  0.9 
o 


K  =  1-0 

o 


R  =  1.1 
o 


1  =  1.2 
o 


i  4 


Fig.  30.  Diagraos  of  twelTe  20-point  inrerse  operators 
designed  for  the  doublet  function  for  R  =  0.25. 
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-0*0000 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0001 

-o.cooo 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0002 

-0.0000 

-0.0000 

-0.0000 

-0.0000 

-0.0001 

-0.0005 

-0.0000 

-0.0000 

-0.0000 

-0.0001 

-0.0004 

-0.0012 

-0.0000 

-0.0000 

-0.0000 

-0.0003 

-0.0010 

-0.0030 

-0.0000 

-0.0000 

-0.0002 

-0.0009 

-0.0029 

-0.0072 

-0.0000 

-0.0002 

-0.0009 

-0.0031 

-0.0080 

-0.0175 

-O.OOOl 

-0.0011 

-0.0039 

-0.0102 

-0.0221 

-0.0424 

-0.0016 

-0.0069 

-0.0169 

-0.0339 

-0.0610 

-0. 1026 

-0.0205 

-0.0437 

-0.0730 

-0.1127 

-0.1686 

-0.2484 

l.OOOO 

1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

0.0797 

0. 1577 

0.2320 

0.3008 

0.3619 

0.4132 

0.0064 

0.0249 

0.0538 

0.0905 

0.1310 

0. 1707 

0.0005 

0.0039 

0.0125 

0.0272 

0.0474 

0.0705 

0.0000 

0.0006 

0.0029 

0.0082 

0.0172 

0.0291 

0.0000 

O.OOOl 

0.0007 

0.0025 

0.0062 

0.0120 

0.0000 

0.0000 

0.0002 

0.0007 

0.0022 

0.0050 

0.0000 

0.0000 

0.0000 

0.0002 

0.0008 

0.0020 

0.0000 

0.0000 

0.0000 

0.0001 

0.0003 

o.oooe 

0.0000 

0.0000 

0.0000 

0.0000 

0.0001 

0.0003 

R  =0.1 

R  *0.2 

K  *0.3 

R  *0.4 

fi  =0.5 

R  *0.6 

o 

o 

o 

o 

o 

o 

-0.0002 

-0.0005 

-0.0010 

-0.0015 

-0.0014 

-0.0012 

-0.0006 

-0.0014 

-0.0025 

-0.0037 

-0.0035 

-0.0030 

-0.0014 

-0.0030 

-0.0053 

-0.0077 

-0.0073 

-0.0064 

-0.0031 

-0.0064 

-0.0108 

-0.0156 

-0.0149 

-0.0132 

-0.0069 

-0.0133 

-0.0218 

-0.0312 

-0.0301 

-0.0273 

-0.0152 

-0.0276 

-0.0440 

-0.0625 

-0.0606 

-0.0561 

-0.0336 

-0. 3575 

-0.0888 

-0.1250 

-0.1222 

-0.1152 

-0.0742 

-0.1198 

-0.1792 

-0.2500 

-0.2463 

-0.2368 

-0. 1639 

-0.2494 

-0.3617 

-0.5000 

-0.4962 

-0.4866 

-0.3618 

-0.5192 

-0.7301 

-1.0000 

-1.0000 

-1.0000 

1.0000 

1.0000 

1.0000 

1.0000 

0.7522 

0.5832 

0.4529 

0.4802 

0.4954 

0.5000 

0.3733 

0.2838 

0.2052 

0.2306 

0.2454 

0.2500 

0.1852 

0.1381 

0.0929 

0. 1108 

0.1216 

0. 1250 

0.0919 

0.0672 

0.0421 

0.0532 

0.0602 

0.0625 

0.0456 

0.0327 

0.0191 

0.0255 

0.0298 

0.0312 

0.0226 

0.0159 

0.C086 

0.0122 

0.0147 

0.0156 

0.0112 

0.0077 

0.0039 

0.0058 

0.0072 

0.0077 

0.0055 

0.0037 

0.0017 

0.0027 

0.0034 

0.0037 

0.0026 

0.0017 

0.0006 

0.0010 

0.0014 

0.0015 

0.0010 

0.0007 

R  *0.7 

R  =0.8 

R  =0.9 

R  =1.0 

R  =1.1 

R  =1.2 

o 

0 

o 

o 

o 

o 

Table  4.  Values  for  twelve  20-point  inverse  operators 
designed  for  R  =  0t23, 
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APPENDIX  C.  Fortran  computer  programs. 


'^WQ  main  programs  and  several  subroutines  were  written  to  perform 
the  routine  computations  necessary  to  accoaplisli  this  thesis.  They  are 
written  for  the  Pennsylvania  State  University  IM  7074  computer. 

This  appendix  briefly  describes  the  purpose  of  each  program  and 
how  they  may  be  used. 


A.  Identification 


Title:  Executive  Program  MINLMAX 
Progranmier:  Kerdler 
Date;  February  1964 

B.  Purpose 

Coordinates  the  functions  of  several  subroutines  that  compute  the 
simplicity  of  a  deconvolved  signal. 

C.  Usage 

1.  Calling  Sequence:  None.  This  is  a  main  program, 

2.  Space  required:  Approximately  4000  locations. 

3.  Input  and  Output  Formats 

a)  Input 

R  noise  to  signal  power  ratio 

DELT  input  signal  sampling  increment 

KUT  total  number  of  different  delay  times  used 

XNL  u -level  increment 

NU  total  number  of  u-levels 

Format  (F10.5,  F5.3,  15,  F5.3,  15) 

number  of  sample  points  in  inverse  operator 
ghost  amplitude  (negative  floating  point  number) 
inverse  filter  lag 

Format  (15,  F5.3,  15) 

LAGS  del*  ime 

Format  (15) 

b)  Output 

Printer 

TEK 

ULEV 
RH 
2 

LAGS 


signal  energy  concentration  for  a  particular 
u-level,  LAGS,  RH  and  R. 
u- level 

ghost  amplitude 
delay  time  in  seconds 
delay  time  in  zeroes 


NF 

RH 

lAGF 


Format  (6X,  E16.8,  5X,  F5.3,  6X,  F7.3,  5X,  F6.3,  5X.  15) 


} 

? 


■» 


50 

100 

150 

?00 

220 

230 

240 

C9<j 


PR5QPAM  M|»j{MAX  r  8  STNthtfTXC  Si3SALS 

input  GrNpBATEs  signal  Sin 

CONVOL  nfGHOSTS  Si  it  fONVOLUTiON  WITH  Pflt 
IiWessp  FILTEP  OF  OOi'PLET  FUNCTION 
mInmAX  determines  BFST  OUTPUT  SIGNAL 
USING  simplicity  CSITFPion 

I^IMFNSi  F  { SO  5  *S  <  SOO  )  *0!  lOQ)  » F!  lOO  )  •T  1 100)  «S#  NO !  500  I  ‘Ci  1000 !  * 
i  BLlMtSOl 

rfad  ?oo.p»oflt»t.ut  *XNL»NU 
CALL  INPUT! s,nS»ph.LAGG) 

OPiNT  100 

REAP  220.NF .PH.LAGP 

IF!NF)909»999.2 

RFAO  250. IF! ! ) ,!=1 .NS) 

00  30  Jal.ifUT 
SUMl-0.0 
00  20  1*1 *NU 
Ef  n»c.o 
T?n*o.o 
RFAO  240. LAGS 

Z.FLOATF!lAGS+n*0FL.T 

CALL  cONVnLIF.S.NF.NS.LAGS.XNL  .NU»F  *T  .U.SIJMI  ) 

00  50  K«1*NU 

IFCTIK 114.4.5 

TFKeFlIf  )/iT»)  )»OFLT*SUMl  ( 

GO  TO  5 
TEK^O.O 
ULEVstTIK) 

PRINT  150.TEF.ULEV.RH.7.LAGS 

FORMAT! IHl .iJX.THMINlMAX.fiX.SMLEVEL.eX.SHG-AMP.EX.SHG-LAG.AX.eHZEP 
1  0ES//1 

FORMAT  !F,X.F16.S»5XtP5,5,AX.F7.5.SX.FA.5,SX»!S  I 
FnRMATIFlO.S^ce.Y.tSjFS.A.IS) 

FORMATns.FS.-^.I*) 

format !4E16. A) 

FORMAT! 15 1 

STOP 


PROGRAM  SlMPTtST  FOR  KARTHi^AKBS 

CONVOL  OEGHOSTS  SU)  PV  CONVOLUTION  WITH  FIIJ 

INVERSE  FILTER  OF  DOUBLET  FUNCTION 

MINMAX  determines  rest  OUTPUT  SIGNAL 

USING  SIMPLICITY  CRITERION 

DIMFNSION  Ff AOI .SIAPOI.UI ICO j  * r I  100  I ♦ T ( l 00 { 

RFAO  ?00*R.0«^LT.rUT.XNL*NU 
CALL  INPUTI5*NS.ph.LAGG> 

PRINT  100 

1  RFAO  ?20*NF.RH»LAGF 
IFINF»999»99R»2 

2  read  230.»FIIJ.I»1»NF> 

DO  30  J-lttf.UT 
SUMl-0.0 

DO  20  I«l»NU 
E(1 1-0. 0 
20  Tl 1 I-O.O 

READ  2A0.LAGS 
Z«FLOATFaAGS+l)*OELT 

CALL  C0NV0LIF.S.NF.NS»LAG5*XNL*NU.F.T,U*SUM1I 

DO  30  X«1 »NU 

IF(T(Xn*»A»3 

3  TEIC»FIX)/ITlK)*OELT*SUMn 
GO  TO  5 

A  Tek-C.O 
5  ULEV«U(K) 

30  PRINT  150tTEK»ULEV.RH,Z*LAGS 
GO  TO  1 

100  FORMAT! IHl tlZX.THMlNlMAXteXfSHLEVFL  fAX ♦5HG“AMP»6X*5H6-LAG*4X*OH?F9 
10ES//I 

150  FORMAT!6XtF16,8#^X,F5,3.AX.F7,3,Ax,F6.3*^X»t5) 

200  format tFlO, 5. F5.3,15»F5. 3.15) 

220  FORMAT! 15  .FA.I.IA  ) 

230  F0RMAT!AF16.fi) 

2A0  FORMAT! 15) 

999  STOP 
END 
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A.  Identification 

Title:  Subroutine  INPUT 

Programmer:  Merdler 
Date;  March  1964 

B.  Purpose 

To  generate  a  synthetic  seismogram. 

C.  Usage 

1.  Calling  Sequence:  Call  INPUT  (S,  NS,  RH,  lAGG) 

Uses  library  functions  SETRND,  RANDF,  and  SAVRND 
stored  in  the  PSU  IBM  7074. 

2.  Arguments  or  parameters 

S  synthetic  seismogram 

NS  number  of  points  in  synthetic  signal 

RH  ghost  amplitude 

LAGG  ghost  delay  time  built  into  synthetic 

3.  Space  required:  Approximately  2500  locations 

4.  Input  and  Output  Formats 

a)  Input 

MS  number  of  points  in  primary  reflectivity 

function 

LAGG  built-in  ghost  delay  time 

LAGB  band-limiting  filter  delay  time  (usually  zero) 

RH  ghost  amplitude  (negative  floating  point  number) 

N  starting  random  number  for  RANDF 

R  noise  to  signal  power  ratio 

Format  (415,  F5.3,  120,  F10.5) 

S(I)  input  primary  reflectivity  function 

BLIM(I)  bancHimiting  filter 

Format  (7F10,5)  for  both 

b)  Output 

S(I)  sjmthetic  seismogram  (4E16.8) 

NS  number  of  points  in  synthetic 

N  final  random  number  generated  (120) 
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D.  Method  or  Algorithm 

This  program  convolves  the  input  priniary  function  S(I)  with 
BLIM(I)  and  measures  tne  total  power  of  S(I)*BLIM(I) .  It  then 
generates  NS  random  numbers  RAND(I)  and  knowing  R  computes 
the  weighting  constant  RN.  Ghosts  of  amplitude  RH'S(I)  are 
added  to  S(I)  with  delay  time  LAGG.  Finally  the  primary  plus 
ghost  synthetic  is  convolved  with  BLIM(I)  and  RN*RAND(I)  is 
added  to  thiu  filtered  signal. 


?t>R  SYNTHETIC  SIGNiLC 


SUBROUT  I Nt  INPUT{S.NS*f»H,LAGGi 
DiMENSTON  f?AHD?500>  jS(500  WG{  lOOC  i  .BLIf-USO) 
READ  100tN*NR,R.NB*LAGR 
CALL  SETRNO(N) 

DO  5  I*1*NR 

5  RANDC I )=PANDE(2,0)-1.0 
CALL  SAVRND(N) 

READ  200*L*G?*LAGG*RH 
NS=NR*?L^<“’R+2) 

DO  10  1*1 *NS 
10  S(I)»0*0 

DO  20  1«1*NR 
K=I+f I-1)*LAGP 
Ls<+LAG6+1 
S(K)*RAN0( 1 1 
20  S(L)«S(L,+RH*RANn(I) 

CALL  SETRND(N) 

DO  25  I=1»NS 

25  RAND?  n*RANDr(2.0}-l,0 
CALL  SAVRNDJN) 

PRINT  300»N 
SUM1«0*0 
SUM2=0.0 
DO  30  1=1 tNS 
SUM1  =  SUM1  +  (S(  m**2 
30  SUM2=SUM2+(RANn( n 
RN=SORTF { R»$UM! /SUM2 ) 

DO  40  1=1 *NS 
40  S{  I  )=5{  n+RN*RANn(  ! ) 

READ  400t  IBLIVK  I  •/ *Iel*NE) 

CALL  C05KlP{RLiM*S*NPtNS*LAGB»G»NG) 

BIGX=0.0 

DO  50  I=1»NG 

IFIAPSFIBIGXJ-ABSECGI I ) 1)45 f 50 #50 
45  BIGX=Gnj 
50  CONTINUE 

DO  60  1=1 tNS 
60  S( 1 .=G( n /ABSFIBIGX) 

PUNCH  500f{Sn)*I  =  l»NS) 

100  FORMAT! I20»15»F10,5*2I5J 
200  F0RMAT(2I5  •F«;,3 
300  FORMAT! 120) 

400  FORMAT!7F10,5) 

500  F0RMATI4E16.F) 

RETURN 

END 


FOH  EARTHQDiUCS  SEISKOQfiAMS 

SUBROUTINE  INPUTCS*NS.RHtLAGG) 
DIMENSION  S(500) 

READ  lOO*NS.LAG6fRH 
READ  ?00* CS( I) .I«l *NS> 

BIGX.0.0 
DO  ID  I-l *NS 

IF(APSF(BIGX)-ABSE{S( 1) M 5* 10. 10 
5  BIGX.Sd) 

10  CONTINUE 
DO  20  I«1.NS 
20  SdleSt  n/ABSFeBIGX) 

PUNCH  500.{SCn.I*l.NS) 

100  FORMAT{2I5.F5.3) 

300  FORMATdSFS.l) 

500  FORMATUE16#8l 
RFTURN 
END 
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A.  Identitication 

Tide:  Subroutliie  COSKIP 


Pfograraner :  Merdler  and  Watson 


Dace.’  January  1964 

B.  Purpose 

To  convolve  two  functions  and  eliminate  multiplying 
zeroes  by  skipping. 


C.  Usage 

1.  Calling  Sequence:  Call  COSKIP (F,S,NF, NS, LAGS. G, KG) 

2.  Arguments  or  Parameters 


F 

first  function  (inverse  operator) 

S 

second  function  (input  signal) 

NF 

number  of  points  in  first 

function 

NS 

number  of  points  in  ; jcond 

function 

LAGS 

spacing  between  points  in 

inverse  operator 

G 

deconvolved  input  signal 

NG 

number  of  points  in  G 

3.  Space  required;  Approximately  2000  locations 

4.  Input  and  Output  Formats 
Determined  by  main  program  MINIMAX 

5.  Method  or  Algorithm 


NG  MIK(L,NF) 

G(I)  ^  Sj  I-(J-l)  - (LAGS +1)1 -FCJ) 

I-i  MAX(1,M)  ‘- 


where  L  l+( I-l) /(LAGS-l) 


M  _  (I-NS)/(iAGS  +  l)  -1 


PROGRAM  COSKIP 
CONVOLVFS  TWO  FUNCTIONS 

fliminatfs  multiplying  ?F?>0FS  py  skipping 

SUBPOUT  INF  COSK!P{F*.S*NF»NS*LAGS*G»NGI 
DIMENSION  FtSO) •SISOO: *Gi 1000) 

NG*^NS+NF+  ( NF-1 )  »L  AGS-l 
00  10  I«1*NG 
L*l+( I-l) /<LA0S+1 ) 

KUP*XMIN0F{L#NF) 

M*t{  I-NS)/(LAGS+1)+1 
B=I I-NSI/(LAGS+1I+1 
iFie-FLOATFIM) )3*3t2 

2  MbM+1 

3  KLOW=XMAXOF(1»M) 

Gn)*0*0  ' 

DO  10  J*KLOW*KUP 
MS=I-( J-1 )*{LAGS+1 ) 

10  G(  I  )*Gn  )+SIMS)*F(  J) 

RETURN 

END 
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A.  Identification 

Title:  Subroutine  CONVOL 

Progransaer:  Merdler 
Date;  March  1964 

B.  Purpose 

CONVOL  convolves  two  functions  using  the  COSKJP  algorithm 
but  does  not  store  G(I). 

C.  Usage 

1.  Calling  Sequence 

Call  CONVOL  (F, S.NF.NS, LAGS,XNL, NU.E, T, U, SUMl) 

2.  Arguments  or  parameters 
F,S,NF,NS,LAGS  same  as  for  COSKIP 
XNL,Nll,U  defined  in  main  program 
EjTjSUMl  are  used  in  subroutine  MINMAX 

3.  Space  requ’ced:  Approximately  600  locations 

4.  Input  and  Output  Formats 
Determined  by  main  program  MINIliAX 

D.  Method  or  Algorithm 

Same  algorithm  used  by  COSKIP 
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A,  Ident t f ica t icn 

Title:  Subroutine  MINMAX 
Prograimaer;  Merdler 
Date;  February  1964 


B.  Purpose 

To  compute  the  normalized  percentage  of  energy  of  tiie  filtered 
signal  G(I)  above  a  specified  u-level  per  unit  of  time. 

C.  Usage 

1.  Calling  Sequence:  Call  MINMAX(G,XNL,NU,E,T,U,SUM1) 

2.  Arguments  or  parameters 

G  filtered  signal  (deconvolved  input  signal) 

E  energy  of  ABSF(F(I))  for  G(I)>u 

T  total  time  width  of  ABSF(F(I))  for  G(I)>  u 

SUMl  total  energy  of  ABSF(G(I))  for  u-level  =  0 

XNL,NU,U  defined  in  main  program  MINIMAX 

3.  Space  required;  Approximately  Vuv.  locations 

4.  Input  and  Output  Formats 
Determined  by  main  program  MINIMAX 

D.  Method  or  Algoritlun 


TEK(I.J) 


T 

^[aBSF(C(I))-U(J)J  ^ 

NG 

T-DELT-  21  CaBSF(G(I))]^ 
I  1 


J  =  1,2, 


,NU 


where  T  =  total  number  of  G(I)'S  for  which  ABSF(G(I)) ^  U(J) . 


J»  O"  l> 


program  convol 

CONVOLVFG  TWO  FUNOTIONF 

ELIMINATES  MULTIPLYING  ZFROES  PV  SKIPPING 
NO  STORAGE  OF  OUTPUT 

SUBROUTINE  CONVOL  <F»StNFtNS»LAGS.XNL»NU»rfT.UfSUMl) 
dimension  FISO) tSISOO) fUl 100) .F( 100 ) .T( 100) 
NG=NS+NF+  <NF-1 )*LAGS-1 
DO  20  1=1 fNG 
L=1+CI-1 )/{LAGS+l ) 

KUPsXMlNOFILtNF) 

I-NS)/(LAGS+1 )+l 
P=CI-NS)/«LAGS+1)+1 

1F(B-FL0ATF(M) 

2  MaM+l 

1  KLOWsXMAXOFd  ♦¥) 

G*0.0 

DO  10  J*XLOW»KUP 
MS=I-C  J-1  )*aAGS+l  ) 

10  G=«G+S(MS)*F(  J) 

CALL  MINMAX{G»XNL»NU»E»T»UtSUMl) 

20  continue 

RETURN 

END 


SUBROt.iTINF  MINMAX(G»XNL»NU*F»T,U»SUM1  ) 

MINMAX  COMPUTES  THE  NORMALIZED  PERCENTAGE  OF  ENERGY 

OF  THE  filtered  SIGNAL  Gd)  APOVE  A 

SPECIFIED  level  U  PER  UNIT  OF  TIME 

DIMENSION  F{50)  tSISOO)  ,UC100)  ♦Ed00)»Td00) 

SUM1=SUM1+G**2 

DO  5  I»lfNU 

U{ I )*XNL*FLOATF{ I ) 

DFLTA»ABSF<G)'U< I ) 

IF{DFLTA)fl,9.6 
Ed  )=E<  I  )+DFLTA**2 
Td)*T{  I)+1.0 
RETURN 
END 


86 


A. 


tden ficat ian 


Title;  Executive  Program  OUTPUT 
Progranaaer ;  Merdler 
Date:  February  1964 
B-  Purpose 

Coordinates  the  functions  of  subroutines  COSKIP  and  PLOTM 
for  visual  display  of  the  correct  inverse  operator,  input 
signal  and  filtered  input. 

C.  Usage 

1.  Calling  Sequence 

None.  This  is  a  main  program 

2.  Space  required:  approximately  3000  locations 

3.  Input  and  ©"atput  Formats 
a)  Input 


NF 

NS 

LAGS 

LAGF 

RH 

DELT 

R 

LAGG 

SCALE 


number  of  points  in  inverse  operator 
number  of  points  in  input  signal 
theoretical  ghost  delay  time 
inverse  operator  lag 
ghost  amplitude 

input  signal  sampling  increment 
noise  to  signal  power  ratio 
actual  ghost  delay  time 
plotting  scale  for  subroutine  PLOTM 


Format  (415,  2F5.3,  FlO.3,  15,  F5.3) 


F(I) 

S(I) 

b)  Output 
Printer 


inverse  operator  4E16.8 

input  signal  4E16.8 


Values  for  NF,R,RH,LAGF 
Values  for  NS, DELT,  LAGG 

Curves  for  F( I)  ,  S(I)  ,  r(I 'i*S(  I)  , RH- (F(  I)*S(I)  ) 
and  F(I)*S(T)-(1-RH) 


Printed  using  PLOTM 


PROGQAM  output  fob  STMTHETIG  SiflHALS 

DIMENSION  FCnOj  ♦SIFOO}  tG?  1000}  ,Hl^OO}  tfl'iOO}  ♦GRAPH f?0) 

1  PPAO  ^0  ♦NF  tNS  *  L AGS  tL  AGF  t  PH*DEL  T  «R  tL  AGfitSCALC 
f F f ) 99^ ^999 , 2 

2  RFAO  1 00  ♦  {  F  <  I )  ♦  I  E  j  ) 

PFAO  lOOt CSJ f ) ♦ I-l tNS} 

CALL  COSif  IP(r,G.NF.N*;,tAGStG*NG) 

DO  10  <■! tNS 

MS*  f  LAGF-i )  #{lags  +  1  1  *r. 

L5«MS*LAGS*l 
THK) ■RH*G{MS) 

10  ElO'GfLSl+HIJt) 

PRINT  300tNFtRtRH.LAGF 
PRINT  310 

CALL  PLOTM{.SF,FtSCAtF» 

PRINT  AOOtNStOELT.LAGG 
PRINT  310 

CALL  PLOTMiNGtStSCALF} 

PRINT  500 
PRINT  310 

CALL  PLOTMlNG.GtSCALF} 

PRINT  600 
PRINT  310 

CALL  PLOTMINStHfSCALE 5 
PRINT  700 
PRINT  310 

CALL  PLOTMINStEtSCALE} 

GO  TO  1 

50  format  (4I5t2F5,3tFlC,TfIAtFA,‘») 

100  FORMATlAEife.gj 

300  format ( IHl tAHNFst  t5t2y,2HRstF7,5t?F.THPM«tF7,A,2Xt5HLAGF« ♦ I A//1h  ♦ 
12?MOPTImUm  InvFRSF  FIlTER/IM  ,4HF n ) ♦ AXt 9hT IMF t sfc ) ) 

310  format  t  IH  ^20*  t*M-l  ,Oti5X  tAH-O,*;  f  17X  t  AHO.Of  1  7X  ♦3MO,5t  1  7X  ♦^m  *0/ 
121X ♦ iHT# ifix t IHT ♦ i9Xt IHTt 19X ♦ ihT t IRX  t IHT ) 

*00  format ( lH0t3HN5=# 15 t2Xt5H0FL7=  tFS.I .1 ^»3hSEC ♦ 1 < ♦ IOhghOST  L  AC», f 5/ 
IIH  tl2HINPUT  SIGNAL/IH  ♦ AHS I  11 » AX tOMT IMF  I  SEC  )  > 

500  Format { iHOt lAHFiLTFpFu  input/ih  fAHon ) tAXtOHiiMEf sfo  i 
600  FORMAT! iHOt 12HGH05TS  ALONE/IH  tAHHI I) ♦ AX ♦ 9HT IMF | SEC ) ) 

700  format nH0t26HFILTFRFn  INPUT  PLUS  GHOSTS/ 

IIH  fAHEI  I  )  ,<,Xt9MTIMF{SFr  )  ) 

999  STOP 
FNO 
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progcam  output  poo  P-QUAKE  data 
MMPNSfON  F«  AOl  .GtADDO)  .GPAPHj  jn) 

I  PFAD  A0.NE«N‘;»tAG5,t  AGP.OM.nrt  T.PtLAGr.tAfALF 
IF{NF)909*999*? 

?  PFAD  IGO.tPfll  «t«l  tN'-l 
ppAo  ICO* (5( n * r«i tNS> 

CALL  rOSKlP<F.5*NF»N5»LAG£.G*NGJ 
PRINT  300.NE«R.RH,LAGE 
PRINT  31C 

CALL  PLOTMUv-.F, scale* 
r>RINT  400*NS.nFLT  .LAGG 
PRINT  310 

CALL  PL0TM(NS.S*SCALF> 

PRINT  500 
PRINT  310 

CALL  PLOTMING.GtSCALF* 

GO  TO  1 

50  format  Ur5t?F5.3*F10.1»fA,F5*-»* 

100  F0RMATUE16*R» 

300  FORMAT? 1  HI  * 3HNF= . 1 5  v ?X .2HR» ,F 7, 5 ♦ 2X » 3mPm. .FT. 3 . 2X .5HLAGF* *  I S/ZIM  ♦ 
122M0PTIM'JH  inverse  F  ILTER/IH  .AHF?  IJ  ,4X*9hTIME(SEC*  I 
310  format? IH  .20X.AH-1,0*15X 5  * 17X * 3HO.O ♦ 1 7X ♦3M0. 5 ♦ 1 7X . 3H1 .0/ 
123X.1HT  tl8XtlHT.19X,lHT.19X»lMTtl9X*lMT ) 

AOO  format? 1M0.3HNS=  *  1 5 . 2X .SmDFL T» * F5, 3  ♦  1 X  * 3HSEC ♦ 1 X ♦ lOHGMOST  LAG* *15/ 
IIM  tl2HINPUT  SIGNAL/IM  .AHS ? I* »AX .RHT IMF ? SEC M 
AOO  FORMAT ? IHOtlAHFILTFRFO  INPUT/IM  * 4HG flJ ♦ aX ♦RMT I Mf ? SFC ) * 

999  STOP 
END 
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A.  Identification 


Title:  Subroutine  PLOTM 


Progrannner :  Merdler 
Date;  January  1964 


B,  Purpose 

To  plot  any  two  dimensional  array.  The  program  normalizes  the 
quantities  involved  prior  to  plotting.  Vertical  scale  is 
chosen  by  user  (width  of  graph). 


C.  Usage 


1.  Calling  Sequence:  Call  PLOTM(K, YORD, S) 

Uses  library  functions  SCAI£F  and  GRAPHF  stored  in  the 
PSU  IBM  7C74. 

2,  Arguments  or  Parameters 

K  number  of  points  to  be  plotted  K<  2000 

YORD  value  of  ordinate 

S  scale  value  (maximum  width  of  curve  <  10  inches) 


How  to  determine  S; 

Suppose  we  desire  that  the  maximum  width  of 
our  graph  to  be  4  inches.  We  know  that  the 
normalized  values  of  YORD  will  range  between 
±1.0. 


Scale 


Max  YORD  -  Min  YORD 
Graph  width 


0.5. 


In  this  case  YORD  will  be  plotted  with  an 
accuracy  of  0.05  units. 


3.  Space  required:  Approximately  2000  locations 

4.  Input  and  Output  Formats 
a)  Input 

Determined  by  main  program 


b)  Output 
Printer 


Vnluev;  of  YORD(I) 

N  mber  of  point  I 
Asterisk  for  value  YORD(I) 

Format  (IHS,  F8.4,  3X,  15,  5X,  20A5) 

5.  Accuracy: 

Values  of  YORD(I)  may  be  plotted  no  closer  than  O.Ol  uni 
Method  or  Algorithm 


Uses 

fPSU 


library  funr  ions  SCALEF  and  GRAPHF 
programming  consultant). 


written  by  Forney 
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\  SUBROUT  INF  PLOTMCK  tYORn»5;) 

2  dimension  graph? 20 ) tVnPDC 2000) 

BIGXsO.O 
4  DO  7  1  =  1  tic 

IF{ARSF(RIGX)~ARSF(YORn{  I  )  )  )6t7t? 

6  BIGX«YORD{I) 

7  CONTINUE 

8  DO  9  1=1 tX 

9  YORDf I )  =  YORD{ I ) /ABSF(PIGX  ) 

10  DO  14  I=ltX 

11  D=SCALEF(-1.0tS) 

1?  DUMMY. GRAPHF {GRAPH tYOPD? I ) #1H») 

13  PRINT  lOOtYQRD? 1) tItGRAPH 

14  CONTINUE 

15  PRINT  200*BIGX 

100  FORMAT! lHStF8*4»5XtI5t5X*20A5) 

200  FORMAT? IH  t21H  NORMALIZING  FACTOR  =tF10.4) 

16  return 

17  END 
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